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ABSTRACT 

Diffusive shock acceleration (DSA) at relativistic shocks is expected to be an important acceleration 
mechanism in a variety of astrophysical objects including extragalactic jets in active galactic nuclei 
and gamma ray bursts. These sources remain good candidate sites for the generation of ultra-high 
energy cosmic rays. In this paper, key predictions of DSA at relativistic shocks that are germane to 
production of relativistic electrons and ions are outlined. The technique employed to identify these 
characteristics is a Monte Carlo simulation of such diffusive acceleration in test-particle, relativistic, 
oblique, magnetohydrodynamic (MHD) shocks. Using a compact prescription for diffusion of charges in 
MHD turbulence, this approach generates particle angular and momentum distributions at any position 
upstream or downstream of the shock. Simulation output is presented for both small angle and large 
angle scattering scenarios, and a variety of shock obliquities including superluminal regimes when the de 
Hoffmann- Teller frame does not exist. The distribution function power-law indices compare favorably 
with results from other techniques. They are found to depend sensitively on the mean magnetic field 
orientation in the shock, and the nature of MHD turbulence that propagates along fields in shock environs. 
An interesting regime of flat spectrum generation is addressed; we provide evidence for it being due to 
shock drift acceleration, a phenomenon well-known in heliospheric shock studies. The impact of these 
theoretical results on blazar science is outlined. Specifically, f ermi-LAT gamma-ray observations of these 
relativistic jet sources are providing significant constraints on important environmental quantities for 
relativistic shocks, namely the field obliquity, the frequency of scattering and the level of field turbulence. 
Subject headings: Cosmic rays: general — particle acceleration — shock waves — 
magnctohydrodynamics — gamma-rays: sources — blazars 



1. INTRODUCTION 

GoUisionlcss magncto-hydrodynamic (MHD) shocks are found 
.io ^diverse environments ranging from the inner heliosphere to 
central regions of distant galaxies and other astrophysical 
kSPjscts. Particle acceleration at these coUisionless shocks is be- 
^'ufeved to be a common phenomenon in space plasmas. In the he- 
^3sphere, direct measurements of accelerated non-thermal ions 
■ ■ and electrons in various energy ranges at the Earth's bow shock 
(e.g. Scholer et al., 1980, Mobius et al., 1987 and Gosling et al., 
1989) and interplanetary shocks (e.g., Sarris & Van Allen, 1974; 
Decker et al. 1981; Tan et al. 1988; Baring et al. 1997) indicate 
energization processes that are intimately connected to shock 
environs. Outside the heliosphere, non-thermal particle distri- 
butions are inferred from observed photon spectra of supernova 
remnants, pulsar wind nebulae, blazars, and gamma-ray bursts 
(e.g. Blandford & Eichler, 1987, and references therein), all 
of which possess supersonic outflows that are readily shocked. 
Commonly, these non-thermal distributions take the form of 
power-law tails that can extend to thousands or millions of times 
the ambient thermal energies of the particles. 

First-order Fermi acceleration, often called diffusive shock ac- 
celeration (DSA), is believed to be the primary acceleration 
mechanism in most coUisionless MHD shocks. This phenomenon 
arises when charged particles interact quasi-elastically with tur- 



bulent fields in the shock layer, and are diffusively transported 
back and forth across the shock, each time achieving a net gain 
in energy on average. Monte Carlo simulations of this process 
(see Jones and Ellison, 1991, and references therein) have had 
great success in modeling shocks inside the heliosphere and com- 
paring them directly with in-situ measurements from various 
spacecraft (e.g. Ellison et al., 1990b; Baring et al., 1997; Sum- 
merlin & Baring, 2006). It is quite likely that this same process 
is responsible for the power-law tails inferred in astrophysical 
shocks, including relativistic MHD discontinuitcs such as those 
believed to be associated with blazars (e.g. see Stecker, Baring 
& Summerlin 2007) and gamma-ray bursts (e.g. see reviews by 
Piran 1999; Meszaros, 2001). 

Early work on relativistic shocks was mostly analytic in the 
test-particle approximation (e.g.. Peacock 1981; Kirk & Schnei- 
der 1987; Heavens & Drury 1988; Kirk & Heavens 1989), 
where the accelerated particles do not contribute significantly 
to the global MHD structure of the shock. Since such sys- 
tems are inherently anisotropic, due to rapid convection of par- 
ticles through and away downstream of the shock, the diffu- 
sion approximation cannot be applied. This renders analytic 
approaches, such as solution of the diffusion-convection Fokker- 
Planck equation, more difficult for ultra-relativistic upstream 
flows, though advances can be made in special cases, such as the 
limit of extremely small angle scattering (e.g. Kirk & Schnei- 
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dcr 1987; Kirk et al. 2000). Accordingly, complementary Monte 
Carlo techniques, first developed for non-relativistic shock appli- 
cations by Ellison, Jones & Eichler (1981), have been employed 
for relativistic shocks by a number of authors, including test- 
particle analyses for steady-state shocks of parallel and oblique 
magnetic fields by Ellison et al. (1990a), Ostrowski (1991), Bed- 
narz & Ostrowski (1998), Baring (1999), Niemiec & Ostrowski 
(2004), Ellison & Double (2004) and Stecker, Baring & Summer- 
lin (2007). It is such a simulational approach that is highlighted 
here; its accessibility to broad dynamic ranges in momenta is ex- 
tremely desirable, providing a niche for Monte Carlo techniques 
in connecting with observations of astronomical objects such as 
gamma-ray bursts (GRBs) and blazars. 

It should be noted that the most comprehensive way to study 
dissipation, acceleration and wave generation in collisionless 
shocks is with particle-in-cell (PIC) simulations, where parti- 
cle motion and field fluctuations are obtained as solutions of 
the Newton-Lorentz and Maxwell's equations. Relativistic PIC 
codes have blossomed to model shocks in applications such as 
GRBs and pulsar wind termination shocks, focusing largely, but 
not exclusively, on perpendicular shocks (e.g. Gallant et al. 
1992; Smolsky & Usov 1996; Silva et al. 2003; Hededal et al. 
2004; Liang & Nishimura 2004; Medvedev et al. 2005; Nishikawa 
et al. 2005; and Spitkovsky 2008). These works have ex- 
plored pair shocks, ion-doped shocks, Poynting flux-dominated 
outflows, and low-field systems with dissipation driven by the 
Weibel instability. PIC simulations are dynamic in nature, and 
rarely achieve a time-asymptotic state. None of these works has 
demonstrated the establishment of an extended power-law that 
is required in modeling emission from gamma-ray bursts and ac- 
tive galactic nuclei, though note the isolated recent suggestion 
(Spitkovsky, 2008; Sironi & Spitkovsky 2011) of a non-thermal 
tail generated by diffusive transport. The general difficulty with 
explicitly seeing acceleration in PIC codes beyond true thermal- 
ization is perhaps due to the severely restricted spatial and tem- 
poral scales of the simulations, imposed by their intensive CPU 
and memory requirements. With the anticipated advances in 
computational capability over the next decade, PIC simulations 
will become a much more powerful tool for probing DSA. For a 
discussion of relativistic shock acceleration, see Baring (2004). 

To date, much simulational work on DSA at relativistic shocks 
has focused on parallel systems (where the magnetic field di- 
rection is parallel to the shock normal) in which particles ex- 
perience frequent small angle scatterings (SAS), as opposed to 
infrequent large angle scatterings (LAS). In the limit of ultra- 
rclativistic shock speeds, for differential particle distributions 
dn/dp = p~'^ , a power-law index of tr ~ 2.23 is realized, as can 
be found analytically (e.g. Kirk et al., 2000) and numerically 
(e.g. Bednarz & Ostrowski, 1998; Baring, 1999; EUison & Dou- 
ble 2004). However, it is not necessary to assume that SAS is the 
dominant scattering mechanism, nor is it warranted in some sit- 
uations: the phase space for the character of small angle scatter- 
ing to be realized shrinks with increasing shock Lorentz factor. 
Moreover, many astrophysical shocks, such as those in blazar 
jets, are either not parallel, or not ultra-relativistic. Clearly, a 
more robust examination of the parameter space is desirable if 
one is to characterize the emission coming from these objects, 
and use it to probe their shocked plasma environments. 

To effect such a goal, here we have extended our Monte Carlo 
DSA code (Summerlin & Baring, 2006) to include shocks of arbi- 
trary speed and obliquity, including the trans-relativistic regime. 
Additionally, we generally presume an electron-positron plasma 
shock, following current thinking on the nature of GRB outfiows 



(e.g. Piran 1999; Meszaros 2001) and blazar jets, though the 
results apply equally well to ion-dominated relativistic shocks. 
The global structure of the shocks is defined via the Rankine- 
Hugoniot relations, solved along the lines of previous expositions 
(e.g. Double et al., 2004). Principal output includes complete 
momentum and angular distributions, at different distances up- 
stream and downstream of the shock. To demonstrate the valid- 
ity of the simulation, and to distinguish its particular character, 
comparisons are made with both theoretical and simulation re- 
sults of other papers (principally Kirk & Heavens, 1989; Kirk et 
al., 2000; Elhson & Double, 2004; Niemiec & Ostrowski 2004). 
More importantly, we expand on these previous works by ex- 
ploring the parameter space for oblique relativistic shocks com- 
prehensively, focusing on the shock obliquity, turbulence levels, 
and parameters encapsulating the microphysics of the turbu- 
lent interactions, as key variables determining the high energy 
power-law index of the particle distribution. 

We find that, in relativistic shocks, unlike in non-relativistic 
shocks, the microphysics of the turbulence becomes an impor- 
tant factor in determining both the value of the power-law index, 
and how many decades in energy particles are accelerated be- 
fore a power law is achieved. Particles undergoing infrequent 
large angle scatterings consistently produce harder power laws 
than their SAS counterparts and take many more decades in en- 
ergy to realize a smooth power-law. It is also apparent that the 
power-law index is critically dependent upon the subluminality, 
versus superluminality, of the shock, as discussed in Sec. U) We 
find that, as do Elhson and Double (2004) and Baring (2004), in 
superluminal shocks, the power-law rapidly becomes softer with 
decreasing levels of turbulence and increasing obliquity, due to 
the difficulty particles have returning to the shock once they 
have crossed to the downstream side. 

In distinct contrast, in the case of subluminal shocks, a de- 
creased amount of turbulence and increased obliquity can ac- 
tually render the acceleration process far more efficient as par- 
ticles undergo the coherent process of shock drift acceleration 
(SDA), where some particles persistently gyrate in the shock 
layer, preferentially gaining energy due to the kinking of the 
magnetic field. In the limit of no cross-field diffusion and a 
de-Hoffmann Teller frame velocity of nearly c, explored theoret- 
ically by Kirk and Heavens (1989) using semi-analytic solutions 
to the diffusion-convection equation, an extremely low value of 
the power-law index around a ~ 1 becomes possible. However, 
with our simulation, we are able to more readily isolate how 
such flat distributions arise. In marginally subluminal shocks 
with SAS operating, a small fraction of high energy particles 
are reflected off the shock by the kink in the magnetic field. For 
those that do reflect, the angular distribution for subsequent 
shock encounters is such that the transmission region is almost 
entirely depleted, resulting in virtually 100% reflection at each 
shock encounter. These particles essentially become trapped 
and are accelerated to very high energies very quickly, before 
they are eventually lost downstream. The extremely low levels 
of turbulence necessary to permit SDA to act unabated almost 
certainly do not occur in Nature, but the effects of SDA can 
be seen to a lesser degree in shocks with more realistic parame- 
ters. In general, it can be concluded that the power-law indices 
in relativistic shocks can sample a broad range, depending on 
the three basic system parameters explored here. After out- 
lining our simulation technique in Sec. [2] and summarizing our 
method for determining the shock jump conditions in Sec. [31 
our results are presented in Sec. HI and then interpreted in the 
context of blazars in Sec. [SJ 
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2. THE MONTE CARLO SIMULATION TECHNIQUE 

The Monte Carlo Simulation technique employed in this pa- 
per closely follows the pioneering work on this method by El- 
lison, Jones & Eichler (1981) and Elhson & Eichler (1984); see 
Jones & Ellison (1991) for a review. It is a test-particle sim- 
ulation that models convection and diffusion of charges in a 
turbulent, shocked flow, complementing the analytic approach 
of Bell (1978) that was extended to the relativistic regime by 
Peacock (1981). It has been successfully applied in a variety of 
environments including, the Earth's bow shock (Ellison, Mobius 
& Paschmann 1990b), interplanetary shocks (e.g Baring et al. 
1997; Summerlin & Baring 2006), the solar wind termination 
shock (see Ellison et al. 1999), supernova remnants (see Baring 
et al. 1999; Baring & Summerlin 2007), and in the regime of 
highly relativistic shocks that is generally encountered in extra- 
galactic contexts (e.g. Ellison et al. 1990a; Ellison & Double 
2004; Stecker, Baring & Summerlin 2007). The code models 
particle gyration about bulk magnetic fields in convecting fluid 
flows, while having their trajectories perturbed by embedded 
hydromagnetic turbulence. The perturbations mediate spatial 
diffusion that permits some small fraction of particles to transit 
the shock front multiple times, kinematically sampling the dif- 
ference in flow speeds on either side of the shock, and thereby 
being accelerated via first-order Fermi (or diffusive) shock ac- 
celeration (see Bell 1978; Jones & Ellison 1991). The code is 
fully relativistic and transitions seamlessly from non-relativistic 
to relativistic flow regimes; it also treats arbitrary orientations 
of the mean magnetic field. 

The simulation space is divided into a distinct number of grid 
zones distributed along the x-axis, which is here defined to be 
the direction normal to the planar shock surface. The bound- 
aries of these grid zones are locations where the bulk properties 
of the fluid (flow speeds, magnetic fields, etc.) can change. The 
values of these fluid properties are specified a priori, and for 
the test-particle implementation of the simulation in this paper, 
have fixed values throughout the simulation runs. In the simu- 
lations presented in this paper a simple step function shock is 
used with only 2 gridzones: one upstream, and one downstream. 
The field and fluid quantities in these two zones are related 
by the fuUy-relativistic, Rankine-Hugoniot jump conditions, as 
discussed in Sec. [3] below. This construction facflitates the gen- 
eralization to non-linear acceleration regimes (e.g. Ellison & 
Eichler 1984; Ellison, Baring & Jones 1996; see also Elhson & 
Double, 2002, for the first treatment of non-linear modification 
of relativistic shocks), where the energetic particles contribute 
to the grid-by-grid specification of MHD quantities constrained 
by energy/momentum fiux conservation. 

Particles are injected isotropically into the system anywhere 
along the x-axis, though usually an upstream injection is 
adopted. The energy distribution of injected particles can be ei- 
ther mono-energetic, a thermal Maxwell-Boltzmann form at any 
temperature (relativistic or non-relativistic) , or a power-law dis- 
tribution in momentum of arbitrary index. For non-relativistic 
shocks with thermal particle injection, the code automatically 
calculates the Rankine-Hugoniot shock jump conditions to as- 
certain the downstream fluid and field vector values. For rela- 
tivistic shocks, the jump condition solution technique is neces- 
sarily more complicated, as described in Sec. 3. This solution 
is accomplished outside the simulation program, and the jump 
conditions are then input manually as initial conditions for the 
simulation runs. The code can also include multiple species of 
charged particles (e.g. treating hydrogenic and pair plasmas. 



and even contributions from helium) besides the test-particles 
in the determination of the jump conditions. After particles are 
injected into the upstream fluid, they are allowed to gyrate in 
the local bulk magnetic fleld, convecting with the fluid, until it 
is determined that a phenomenological scattering occurs. 

The effects of magnetic turbulence are simulated by specifying 
a local fluid frame mean free path for particle diffusion, given 
by 

A = Ao ( — ) oc p" , Ao = 777-31 , (1) 

where rg = pc/{qB) is the gyroradius of an ion or electron of 
momentum p — mv , mass m , and charge g in a magnetic field 
B = |B| . Also Tgi = muixc/{qB) is the gyroradius of an ion 
with a speed v equal to the velocity component, uix , of the far 
upstream flow normal to the shock plane; here x denotes the 
direction normal to the shock. Without loss of generality, the 
mean free path scale Ao is set proportional to rgi with constant 
of proportionality rj defined via Eq. ([T} . This phenomenologi- 
cal prescription for scattering was adopted in numerous papers 
outlining results from the Monte Carlo technique, including Elli- 
son, Jones & Eichler (1981), Ellison, Jones & Reynolds (1990a), 
Ellison, Baring & Jones (1995, 1996), and Stecker, Baring & 
Summerlin (2007). Following this and other previous Monte 
Carlo work, for simplicity, we set a — 1 , a specialization that is 
appropriate for traveling interplanetary shocks: see Ellison et al. 
(1990a,b), Mason et al. (1983) and Giacalone et al. (1992) for 
discussions about the micro-physical expectations for a . The 
simulation can easily accommodate other values of a , however, 
the spectral results are somewhat insensitive to the choice of 
this parameter — its dominant effect is to modify the relative 
scale lengths for diffusion at different particle momenta. Since 
X > Tg is required for physically meaningful diffusion resulting 
from gyro-resonant wave-particle interactions, the a = 1 case 
is also motivated on fundamental grounds. The mean free path 
represents the spatial scale in the local fluid frame on which 
the momentum vector is deflected by 7r/2, on average. Note 
that for diffusion that is driven by non-gyroresonant interac- 
tions with field turbulence, perhaps grown via filamentation or 
Weibel instabilities, it is quite possible to sample rj < I regimes, 
especially when the ambient magnetic field is quite low. Diffu- 
sion in this domain resembles the Bohm limit of 77 = 1 for gy- 
roresonant diffusion, and accordingly the distributions for shock- 
accelerated charges are only mildly dependent on 77 when it is 
less than unity. For high Alfvenic Mach number shocks, the 
scattering is approximately elastic in the fluid frame, i.e., |p| 
is conserved in this frame for interactions with fleld turbulence 
that perturb a particle's pitch angle 9 , gyrophase and orbital 
gyroccnter. 

When the Alfvenic Mach number Ma is low, the Alfvcn 
waves move with appreciable speed in the fluid frame, so that 
partial inelasticity in scatterings arises. This yields second or- 
der, stochastic diffusion contributions. While these can be rou- 
tinely modeled in the simulation, inspection of Eqs. (8) and (10) 
of Pryadkho & Petrosian's (1997) quasi-linear stochastic accel- 
eration formalism clearly indicates that the stochastic contri- 
bution to the spatial diffusion coefficients is smaller than the 
first order Fermi one by the order of 1/J^\. For the effi- 
cient generation of the high energy power-law tails that are 
the primary focus of this paper, the astrophysical shocks of in- 
terest generally have large enough Alfvenic Mach numbers to 
neglect the effects of second-order acceleration. However, note 
that for near-luminal shocks at slightly suprathermal energies. 
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particles are generally unable to convect upstream against the 
downstream flow and are inexorably swept downstream. In this 
energy regime, other mechanisms acting in the shock environs 
such as second-order Fermi acceleration or electrostatic cross- 
shock potentials may noticeably broaden/heat the downstream 
distribution function. This can then enhance injection into the 
first-order acceleration process, and thereby affect the normal- 
ization of the power-law tail that results, particularly in cases 
of strongly-inhibited injection. Notwithstanding, treatment of 
stochastic acceleration effects will be deferred to future work. 

The simplest invocation of scattering is to isotropize the fluid 
frame momentum over the surface of the sphere in momen- 
tum space (Ellison ct al., 1990a). This is large angle scatter- 
ing (LAS), and physically corresponds to large magnetic dis- 
turbances that completely disrupt trajectories of particles. To 
model moderate or even smaller disturbances, each scattering 
event can be restricted to a much smaller solid angle, i.e. can 
be isotropic on a conical sector of a momentum sphere. The 
angular extent of this spherical sector SOmax becomes an ad- 
ditional parameter for the diffusion. Then multiple scattering 
events are required to realize a full mean free path. This is 
the scattering construct that is employed in this paper. The 
relationship between 56max and A was originally developed in 
Ellison et al. (1990a), but is more succinctly presented in Ellison 
and Double (2004) via 

S9max = , (2) 

where rg is the gyro-radius, and N is the number of times 
per gyroradius the particle is scattered. The limit of small an- 
gle scattering (SAS) corresponds to TV ^ 1 , for which the in- 
crement dp in momentum in a scattering satisfies |(5p|/|p| ~ 
SO max ■ In practice, as will become evident below, for relativis- 
tic shocks the SAS domain is realized when the scattering angle 
satisfies S9max ^ 1/ri , where Fi = l/^/l — u\^/c^ is the bulk 
Lorentz factor of the upstream fluid in the shock rest frame. 

Cross-field diffusion emerges naturally from this scattering 
mechanism since, at every scattering, the particle's momentum 
vector is shifted in the local fluid frame, with the resulting ef- 
fect that the gyroccnter of the particle is shifted randomly by 
a distance of order rg sin 6 in the plane orthogonal to the lo- 
cal field. Transport perpendicular to the field is then governed 
by a kinetic theory description, so that the ratio of the spatial 
diffusion coefficients parallel (k|| = Xv/Z) and perpendicular 
( ) to the magnetic field is given by k±/k\\ = 1/(1 + rf') 
(see Forman, Jokipii & Owens 1974; Ellison, Baring & Jones 
1995, for detailed expositions). Hence, ry couples directly to 
the amount of cross-field diffusion and is a measure of the level 
of turbulence in the system, i.e., is an indicator of (SB/B) . The 
quasi-isotropic diffusion case of 77 = 1 constitutes the Bohm dif- 
fusion limit, presumably corresponding to (SB/B) w 1 . 

As will become clear in Sees. 14.2 1 and 14. 3[ in oblique relativis- 
tic shocks, the resulting energy spectrum is critically dependent 
upon both r] , due to the necessity of cross-field diffusion, and 
the scattering angle SOmax , due to beaming effects, producing 
a broad range of power law indices. For small angle scattering 
(SAS) regimes, SOmax < 1/ri , there is little variation in the 
power law tails when other parameters are held constant, since 
the scatter angle is now less than the relativistic beaming angle, 
and the diffusion process becomes insensitive to the scattering 
kernel. Except for Sec. 14.41 SAS is deployed throughout this pa- 
per. Examples of the differences between small and large angle 



scattering in relativistic shocks can be seen Fig. 2 of Stecker, 
Baring & Summerlin (2007), and also in Fig. [T^ 

In between each of the N scatterings per mean free path, 
the code calculates shock frame gyro-orbit trajectories using 
a semi-analytic solver rather than the more popular Bulirsch- 
Stoer method (Stoer & Bulirsch, 1980). Using the properties of 
the magnetized fluid, the shock frame position as a function of 
time is easily derived analytically. The particle is then moved 
along this analytic trajectory until one of two conditions is met: 
(a) the particle scatters or (b) the particle reaches the edge of 
a grid zone. The solution for the time it takes a particle to 
reach the edge of a grid zone must be performed numerically, 
since it involves roots of a transcendental equation of motion in 
the shock frame - the simulation employs a standard bisection 
technique for this purpose. When a particle crosses a grid zone 
boundary, the local fluid properties change, and the trajectory is 
recalculated and the propagation continues. When distances be- 
tween scatterings are many gyro-radii, the semi-analytic method 
can go from one scattering to the next in one step covering many 
gyro-orbits in a single computational step. The Bulirsch-Stoer 
method will always require at least several steps per gyro-orbit 
due to the curvature of the trajectory. However, if particles scat- 
ter many times in one gyro-radius, the increased overhead of the 
semi-analytic method makes it slower than the Bulirsch-Stoer 
method, but not unreasonably so. 

Particles that do not immediately return to the shock may 
isotropize in the downstream reference frame once they have 
traveled, on average, one mean free path. At this juncture, 
an analytical formula developed originally by Bell (1978) and 
later shown to be applicable to relativistic shocks by Peacock 
(1981; see also Jones & Ellison 1991) can be used to calculate 
the probability Pr that a particle heading downstream through 
a, y- z plane at a particular distance x downstream will return 
upstream of this plane: 

In the above equation, u is the local downstream flow speed, 
and Vf is the speed of the particle in this fluid frame. Particles 
that are deemed to fail to return are removed from the system. 
For those ascertained to be returning, their vector velocity com- 
ponents are also determined probabilistically. The particles are 
isotropic in the local fluid frame and have constant energy in 
the downstream frame of reference thanks to the elastic scatter- 
ing off magnetic turbulence. So, the probability of a particle of 
a given fluid frame momentum returning with a particular an- 
gle cosine with respect to the shock normal, /is , can be found 
for arbitrary values of the particle speed and downstream flow 
speed. The details of this calculation and final result can be 
found in the Appendix, specifically Eq. (|A7|) . Employing this 
result, a simple acccpt-rcject method (Garcia, 2000, Ch. 11) 
can be used to select a value for /ig for particles determined 
to have returned. This statistical decision algorithm circum- 
vents excessive computations of extensive downstream diffusion 
that are irrelevant to the acceleration process; accordingly, it 
speeds up the simulation dramatically. Using the correct angu- 
lar distribution of returning particles, i.e. Eq. (jA7p . is essential, 
guaranteeing that the complete distribution function of particles 
anywhere upstream of the probability of return plane is inde- 
pendent of the choice of x , provided a; > A and isotropy in the 
fiuid frame is satisfied at x . 
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For simulation output, accounting of particles in distinct mo- 
mentum bins is documented. As a result of statistical losses 
in the downstream region, when less than half of the particles 
originally recorded in a given momentum bin are retained, the 
remaining particles are "split" into two particles each with half 
the "counting" value of the original. This technique of parti- 
cle splitting allows the simulation to maintain good statistics 
over a large energy range. This extensive energy range is one 
of the primary advantages that the Monte-Carlo technique has 
over other types of simulations. Compared to hybrid plasma 
simulations and particle-in-cell (PIC) simulations, Monte-Carlo 
simulations are computationally inexpensive, allowing the sim- 
ulation to be run long enough for particles to be accelerated to 
very high energies, well above that of the incoming upstream 
ballistic flow, in a reasonable amount of time. 

In the test particle implementation employed here, the char- 
acteristics of the shock and the functional form of the turbu- 
lence are specified a priori. The test particle approximation is 
entirely appropriate unless there is a significant fraction of the 
total energy present in the accelerated particles. Since the dis- 
tribution of these particles is only known after the simulation 
accelerates them, the shape of the shock can not be adjusted 
to account for their existence until after the simulation is run. 
For non-relativistic shocks, Ellison and Eichler (1984) developed 
a feedback loop technique where subsequent runs calculate the 
modified hydrodynamic shock structure, based upon the accel- 
erated particle distributions of the previous iterations; this non- 
linear acceleration method is not employed here. Also, since 
the choice of the scattering mechanism can affect both injection 
and acceleration of particles, it can strongly impact non-linear 
modifications for relativistic shocks. The influence of different 
scattering scenarios in such non-linear acceleration systems will 
be the subject of future work. 

This implementation also does not retain accounting of the 
amount of time the particle would have spent downstream of 
this "return" plane. In the event that acceleration time informa- 
tion is needed, a retrodictive approach described first in Jones 
(1978) and later applied directly to a Monte-Carlo simulation 
in Ellison et al. (1990a) can be used. One important finding 
is that the interplay between energy boosting and time dilation 
effects leads only to modest changes (Baring 2002) in the accel- 
eration time at plane-parallel relativistic shocks compared with 
standard non-relativistic shock formalism (Forman, Jokipii & 
Owens 1974). The consideration of particle acceleration times 
is beyond the scope of the present work, and will be deferred to 
a future investigation. 

3. MAGNETOHYDRODYNAMIC JUMP CONDITIONS FOR 
OBLIQUE RELATIVISTIC SHOCKS 

In the case of relativistic shocks, the shock jump conditions 
are considerably more difficult to solve than the non-relativistic 
solutions presented in Decker (1988), due to the impact of length 
contraction and time dilation effects on the structure of the six 
conservation equations. There are different approaches to solv- 
ing the Rankine-Hugoniot jump conditions in relativistic MHD 
discontinuities, surveyed in Double et al. (2004; see also Gerbig 
& Schlickeiser 2011, for a recent exposition). Our approach here 
builds upon previous work by Ballard & Heavens (1991) that for- 
mulates the Rankine-Hugoniot conditions in the de Hoffmann- 
Teller frame (de Hoffmann and Teller, 1950; hereafter HT) in a 
manageable form. The HT frame is a shock rest frame in which 
there are no u x B drift electric fields. This can be obtained 
from the local fluid frame by boosting along B, but can also be 



obtained as a combination of two boosts along the axes of the 
coordinate system to avoid a rotation of the coordinate system. 
The system of equations is then transformed from the HT frame 
into the normal incidence frame (NIF, in which the upstream 
plasma flow is parallel to the shock normal or x direction), ar- 
riving at a system of three comparatively simple simultaneous 
equations in which the terms that become imaginary in super- 
luminal shocks are no longer present. These three equations are 
solved numerically after the Jiittner-Synge equation of state is 
invoked to connect key thermodynamic quantities such as pres- 
sure and enthalpy, to the temperatures of the upstream and 
downstream relativistic Maxwell-Boltzmann distributions. This 
method encompasses a broad range of shock conditions, specif- 
ically ranges of sonic and Alfvenic Mach numbers, and transi- 
tions seamlessly from subluminal to superluminal regimes. Our 
results are compared directly with that of the work by Double et 
al. (2004), highlighting similarities, and also differences that re- 
sult from a specific approximation to the downstream equation 
of state employed in that work. 

Before embarking upon the construction and reduction of the 
jump conditions, a brief summary of the subscript conventions 
adopted here for the different frames of reference is given. The 
'f ' subscript will denote a quantity defined in the rest frame of 
the upstream (subscript 1) or downstream (subscript 2) plasma. 
HT frame variables will be subscripted with an 'HT.' To distin- 
guish NIF frame quantities from those measured in the fluid 
or HT frames, they will denoted by an 'S' subscript for the 
shock frame. Additionally, Ob will always refer to an angle the 
magnetic fleld B makes with the shock normal, and 9^ will 
refer to the angle a plasma flow makes with the shock normal. 
When the HT frame is found via a single boost along the di- 
rection of the magnetic fleld, the fleld components are identical 
in the local fluid and HT frames, often the 'f and 'HT' sub- 
scripts will be explicitly omitted for compactness of notation, 
i.e., Bi = Bit = BiHT , etc. 

The character of the solutions to this system of equations is 
controlled by two key parameters, basically the relative scaling 
of the upstream thermal energy or pressure Pi , and the fluid 
frame magnetic fleld energy density _B^j/(87r) to the upstream 
ram pressure piu\^ . Here, uix ~ /5ixsc is the velocity com- 
ponent of the upstream fluid normal to the shock, in the NIF. 
Accordingly, we deflne these via the sonic (A^s) and Alfvenic 
[Ma) Mach numbers: 

These are conventional definitions for non-relativistic shocks, 
and their extension to oblique discontinuities and relativistic 
systems does not lead to unique choices. For example, subjec- 
tivity is involved in deciding between ui^ and ui , and similarly 
for Bixt versus Bu . Hence, we adopt the above definitions (as 
did Double et al. 2004), for which 7^1 is the upstream adiabatic 
gas index, discussed further in Sec. 13. 3[ so that 7giPi/pi is the 
square of the upstream sound speed. 

3.1. The de Hoffmann- Teller Frame Solution 

For subluminal fiows, where Pix/ cos ©en < 1 , the HT frame 
is an obvious choice in which the shock jump conditions can be 
written, since therein the jump conditions reduce to a simple 
form because the fiuid fiows along the magnetic field lines and 
there is no u x B electric field. For the time being, we will re- 
strict considerations to these types of shocks and later trivially 
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Normal Incidence Frame de Hoffmann-Teller Frame 




X 



Fig. 1. — The geometry in tiie normal incidence (NIP; left panel) and de-Hoffmann Teller (FIT; right panel) shock rest frames. Upstream flow speeds 
in the two reference frames are related by Wiht = /^ihtc = Mi/cos©Bfi ■ Upstream and downstream quantities are denoted by subscripts f and 2, 
respectively. In general, the NIF field angle ©bsI differs from the fluid framc/FIT frame value ©Bfi = ©bhti i and likewise for the downstream angles. 
Also, in cases where the FIT frame is obtained by boosting along the field fluid frame direction, the shock plane in the HT frame is rotated from that in 
the NIF due to relativistic aberration effects. For the two-step fluid-to-HT frame transformation protocol adopted here, the shock planes in the NIF and 
HT frames are coincident, i.e. no rotation is involved. 



generalize the results to include superluminal shocks. Four of 
the shock jump equations are defined by the conservation of the 
mass, momentum (2 components), and energy fluxes across the 
shock interface are conserved. The remaining two derive from 
the electromagnetic field constraints V-B = and VxE = , 
the latter of which is trivial in the HT frame, because ^ = 
identically everywhere. 

The form of these jump conditions in the HT frame has been 
derived previously by Ballard and Heavens (1991). Those equa- 
tions are reproduced below with the notable exception that the 
subscript 'y' used in their paper has been replaced with the 
subscript 'z' to avoid confusion when comparing to other works. 
Here, the a; -direction defines the normal to the shock plane in 
the HT frame, the magnetic field everywhere lies in the x - z 
plane, and the y-axis defines the direction oi u x B drift ve- 
locities. All quantities save F , /3 , and B are defined in the 
frame where the plasma is stationary hereafter referred to as 
the "local fluid frame" or "upstream/downstream rest frame". 
For the present, those 3 quantities are defined in the HT frame. 
Setting c = 1 , as is done throughout this paper, conservation 
of mass or particle number fiux along the shock normal gives 

^iPlxPl = r2/32a:P2 (5) 

where pi denotes mass density, and subscripts 1 and 2 denote 
upstream and downstream quantities (labelled by i in general). 
Throughout this subsection, HT subscripts will be omitted, but 
implied. Also, /3 is the flow speed written as a fraction of the 
speed of light, and F = 1/ -^/l — — /3f is the Lorentz factor 
associated with the flow speed /3 . Conservation of the x and 
z components of momentum fiux gives, respectively. 



Ff wi + Pi 



BixBi 
An 



B2xB2z 

47r 



(6) 



This corrects an obvious typographical error in Eq. (26) of Bal- 
lard and Heavens (1991) in their terms involving the enthalpies 
Wi = ei + Pi . The internal energy , which includes the rest 
mass energy density, can be related to Pi and pi through an 
equation of state, as is addressed in Sec. 13.31 below. In the HT 
frame, conservation of energy flux is simply 



T\(i2x W2 



(7) 



Here, the magnetic field contributions to the stress-energy ten- 
sor (see, for example, Eq. (21) of Double et al. 2004) cancel to 
zero precisely because of the pair of equations 



Plx 



Biz _ 
Bix 



tan Ot: 



P2x 



^ = tanGi: 



(8) 

that define the specific choice of the HT frame. The absence 
of such magnetic terms in the energy fiux, combined with the 
compact nature of the momentum fiux conditions, underline the 
attractive simplicity of adopting the HT frame (compare, for ex- 
ample, with the electromagnetic stress tensor contributions to 
the momentum fluxes in Eqs. (25) and (26) of Double et al. 
2004). The trivial VxE ~ can be efiminated, efii'ectively 
being replaced by the HT frame definitions in Equation Fi- 
nally, the Maxwell equation \7 B = defining the absence of 
magnetic monopoles gives 



Bix 



Bo 



(9) 



unaltered by relativistic generalization because it is intrinsically 
relativistic. 

Following Ballard and Heavens, Eqs. ((8)) and ([9|) can be used 
to eliminate the z-components Biz and B2z , and the down- 
stream x-component B2x ■ Their solutions were defined in 
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terms of two upstream parameters A = Y\vj\l p\ and C = 
Bix/i'^^PilgiPix] ■ Here, as an alternative listing, we observe 
that the ratio C/A appears repeatedly in the resulting subset 
of processed equations, so we define this ratio via 



cos^ Qv 



(10) 



which, as a relativistically-modified ratio of magnetic to ther- 
mal (plus rest mass) energy density, is essentially an adap- 
tation of the inverse of the upstream plasma beta parameter 
(3p = SnPi/Bf to oblique, relativistic MHD flows. The second 
prescription for ■0 uses the Alfvcnic Mach number definition in 
Eq. (|4]), together with identity of total magnetic fields in the 
fluid and HT frames, i.e. Siht = . The energy flux equa- 
tion is most easily manipulated, dividing Eq. ([7|) by the mass 
conservation in Eq. dS]): 



Pi 



T2W2 

P2 



W2 



P2\ l~PL-P, 



(11) 



This is just the constant A employed by Ballard & Heavens 
(1991). Next, dividing the z-component of momentum conser- 
vation in Eq. © by Eq. ([7]) solves for /322 : 



/3: 



2z 



P2x - -0 



(12) 



This can be inserted into Eq. (fTT]) . eliminating /32z . Observe 
that viable jump conditions arc only realizable when -0 < I32x ■ 
This is equivalent to requiring that the total Mach number be 
greater than unity. Finally, the x-component of momentum con- 
servation in Eq. ^ can be divided by Eq. ([7]), producing 



Plx + 



hx + 



Pi 




Plz 






\Plx 


P2 




(Plz 


ri/3la;Wl 




\Plx 



(13) 



Here, expressing the ratio Piz/Pix = tanOenTi in terms of 
the de Hoffmann- Teller field angle Obhti , a constant for the 
shock structure, yields an alternative algebraic form. Observe 
also that the second term on the right hand side is propor- 
tional to P2/W2 multiplied by u'2/wi ; the second factor can 
be expressed using Eq. (fTTj) . and the first is a function of the 
downstream temperature T2 through the equation of state, to 
be defined in Sec. 13.31 

Eq. pT|) . with Eq. inserted, and constitute a sys- 
tem of two simultaneous equations with unknowns P2 , W2 and 
/32a; • However, W2 will be related in section 13.31 to P2 via an 
equation of state, rendering the system numerically solvable. 
This set of equations is simple and elegant, being virtually as 
compact as the system for jump conditions at relativistic, plane- 
parallel, hydrodynamic shocks (e.g. Blandford & McKee 1976). 
However, their validity is technically restricted to subluminal 
regimes where the HT frame formally exists. Therefore, to re- 
alize broader applicability, it is necessary to transform them to 
the normal incidence shock rest frame, and thereafter explore 
their numerical solution. 



3.2. Transforming to the Normal Incidence Shock Frame 

In subluminal cases where the HT frame exists, one can de- 
fine a boost velocity j3t in the z direction that transforms from 
the normal incidence frame into the HT frame. The two key 
input quantities in this regard are /3ixs , the shock speed in the 
upstream fluid frame, and Oan , the angle between the shock ve- 
locity and the magnetic field vector in the upstream fluid frame. 
A third parameter that is derivative of these two is the HT frame 
field angle Obhti • As discussed by Kirk & Heavens (1989), 
there is a lack of uniqueness in defining field and flow angles in 
the de Hoffmann- Teller frame, up to rotations. Here, we will 
adopt the following sequence of boosts to effect Lorentz trans- 
formation to the HT frame from the local fluid frame: this will 
be performed by first boosting by /3xsi along the shock normal 
to the NIF, and then boosting by Ptz in the shock plane to ar- 
rive at the HT frame. The planes of the shock in both the NIF 
and HT frames are thereby coincident. This yields a convenient 
definition of Obhti (and Obht2 ) , and is the preferred protocol 
for our simulation due to the enhanced simplicity it permits for 
modeling particle convection and gyration in the shock layer. 
However, it should be emphasized that a single boost along the 
field vector B from the fluid to HT frames yields an aberration 
of the shock plane: it is rotated relative the NIF shock plane, 
as described in Ballard & Heavens (1991), and is illustrated in 
Fig. [TJ Such a rotation leads to a need to account for it when 
deflning fleld and flow angles with respect to the shock plane, 
an unnecessary complication. The two-step fluid-to-HT frame 
transformation approach adopted here was also the preference 
of Kirk & Heavens (1989). 

The flow velocities in the NIF and HT frames of reference are 
related via standard Lorentz transformations 



rt(l+A/3.s) 
(l + /3t/3zs) ■ 



(14) 



These relations can be applied to both the upstream and down- 
stream sides of the shock; their subscripts 1,2 have been sup- 
pressed here for the sake of compactness. In the upstream re- 
gion, where the deflnition of the NIF requires that /3zs = , 
these equations distill down to /3xht = /^xs/r* and /3zht = Pt , 
respectively. Subsequently, taking the ratio of these two up- 
stream equations, one can express the boost speed Pt and 
Lorentz factor Ff = (1 — {3^)^^/^ in terms of /3ixs and Obh 



/3ixs tan Ob HTl 

Pixs tan Obhti 
1 + /S^^s tan^ Obhti 



(15) 



Since, for flux conserving jump conditions in MHD disconti- 
nuities, the HT frame is identical for both upstream and down- 
stream locations, it can be inferred that the identities in Eq. (jl5p 
can be written also in terms of downstream quantities, merely 
via the substitutions /3ixs /32xs and Obhti — >■ Obht2 • 

The relationship between the magnetic field components in 
the two frames of reference is similarly routinely derived via 
standard transformation equations for electromagnetic fields: 



^xS 



B, 



(16) 
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noting that the equation for -Bxht is only one part of the fuU 
transformation equations, which also transform the NIF drift 
electric field in the y direction exactly to zero in the HT frame. 
The ratios of the equations in Eq. (jl6p then simply yield 



tanOBHTi = TttanOBsi 



tan 0t: 



= r, tan Op 



(17) 



for the upstream and downstream HT frame field angles to the 
shock normal. These are recognizable as aberration formulae 
for the electromagnetic/photon field, with the NIF frame field 
obliquity always being less than that in the de Hoffmann- Teller 
frame. Combining this result with the Eq. (fTS)) yields the rela- 
tionship 



/3t = /3ixstan9Bsi 



(18) 



removing a reference to the HT frame from Eq. ([T4|) . The sub- 
luminal condition for the existence of the HT frame, written in 
terms of NIF quantities, is then /3ixs tan 0bsi < 1 • 

To close this sequence of boost algebra, one needs the relation 
between field angles in the fluid frames and the NIF. These are 
derived in the same manner as Eq. ([TT]). yielding 



tan Qv. 



where Fi 



Fin tanOj 



tan Or 



Fos tanOr 



(19) 



1/Vl - Axs : and Fa 



- PLs - Pis ■ 

As a result, Eq. (jl8p can be rewritten using tan0BHTi = 
FisFf tan Qbh to yield a boost speed expressed entirely in terms 
of input quantities: 



FisAxs tanOE 



(20) 



This then routinely rearranges so that the subluminal /3t < 1 
condition becomes the familiar /3ixs / cos Qmi < 1 • 

By replacing the HT frame quantities in Eqs. PT|) . ([T^ and 
(IT5|) with their shock frame equivalents via Eqs. ((HI), (|17p . and 
(IT5|) . the Rankine-Hugoniot relations then become a system of 3 
equations with unknowns P2 , 11)2 , /32xs and /32zs that possesses 
a non-singular mathematical character in the NIF frame at the 
luminal interface /3ixs / cos Qbh = 1 • The system now retains 
only information about upstream and downstream fluid frame 
field angles and thermodynamic quantities, and the transforma- 
tion velocities required to get to the NIF frame from the fluid 
frames. It must be emphasized that an attractive characteristic 
of this methodology is that significant cancellations remove any 
terms that become imaginary or unphysical in a supcrluminal 
shock, revealing a smooth mathematical transition of solutions 
from subluminal to supcrluminal regimes. The specification of 
an equation of state that relates P2 to W2 closes this system, 
rendering it amenable to numerical solution. 

3.3. The Equation of State 

Assuming there are no shear stresses and axial symmetry 
about the magnetic field, the pressure tensor is diagonal. One 
can then form an isotropic pressure P = (P|| -I- 2P_l)/3 , where 
P|l and P± denote the pressure components respectively par- 
allel to and perpendicular to the mean magnetic field. Then, 
the "adiabatic" gas index 7^ , the approximate ratio of specific 
heats, can parametrize the equation of state via the adiabatic 
expansion law for an ideal gas: 



P\/T9 = constant. 



(21) 



Here 7^ ranges between 5/3 for a non-relativistic, compressible 
gas, to 4/3 for an ultra-relativistic gas. With the specification 



of this index, the internal (thermal) energy density P/ (7^ — 1) is 
related to the pressure, so that the total internal energy density 
plus the rest mass energy density can be written 



P 



la 



1 



+ P 



(22) 



where, again, c=l has been assumed, as will be done through- 
out the rest of this work. Reintroducing the subscripts i = 1,2 
to label upstream and downstream fiuid frames, this leads to 
the forms for the enthalpies that are deployed in the Rankine- 
Hugoniot relations: 



P. 



IgiPi 



1 



(23) 



The particular values of the 7^^ can be expressed as a moment 
of the fiuid frame particle momentum distributions upstream 
and downstream, and so can apply to both thermal and non- 
thermal populations. While they are simply prescribed in the 
non-relativistic and ultra-relativistic asymptotic cases, a more 
precise formulation is required to treat the mildly-relativistic 
domain. 

Here it is assumed that the background plasma possesses a rel- 
ativistic thermal Maxwell-Boltzmann distribution that defines 
the Jiittncr-Synge equation of state (e.g. Synge 1957). Then 
the temperature T can be the sole thermodynamic parame- 
ter, and all other thermodynamic quantities can be prescribed 
in terms of it. The equation of state depends on the number 
of species, their masses, and the state of thermal equilibrium 
between the various species, i.e the temperature equilibration 
or otherwise. For simplicity, here a single component plasma 
is adopted, appropriate for an electron-positron pair plasma as 
might be encountered in relativistic jets in extragalactic sources 
such as gamma-ray bursts or blazars. Equations of states for 
electron-ion and other mixed species gases are addressed in Bal- 
lard & Heavens (1991). For a pair plasma, the enthalpy can be 
written in terms of modified Bessel functions: 



Pi 



Rin) + n , R{t) 



3r 



Ki{1/t) 
K2{l/r) ' 



(24) 



where the Ki s are modified Bessel functions of the second kind 
(e.g., see pp. 708-715 of Arfken and Weber, 2001), and 



El 



i = 1,2 



(25) 



is the dimensionless pair temperature (in units of c = 1 ) . This 
equation of state can treat arbitrary sonic Mach numbers, in 
distinct contrast to the approximation employed by Double et 
al. (2004), discussed below, that uses kinematics in high A4s 
shocks to specify the downstream pressure. 

A modest disadvantage of the Jiittner-Synge equation of state 
lies in the complexity of the Bessel function; it is not conducive 
to either analytical or numerical solutions of Eqs. pT|) . (fT2)) 
and ([T3| . However, noting the asymptotic behavior of Eq. ([24| . 
namely R{t) ^ 3t as t — >■ 00 , and R{t) — > 1 + 3t/2 as 
T — >■ , a remarkably good approximation for the function R{t) 
is given by a Fade approximation of 3rd-order: 



i?(r) 



2 + 7t + 12t^ 



6t^ 



2 + 4t + 2r^ 



(26) 
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Fig. 2. — Rankine-Hugoniot relations for r , the compression ratio, 0Bf2 , the angle the downstream fluid frame magnetic field makes with 
the shock normal, and Ou2 , the angle the downstream flow makes with the shock normal, as functions of the upstream NfF rapidity of the 
shock, Fi/Ji = Fis/SixS • Solutions are displayed for various values of the Alfvenic {A4a) and sonic {Ms) Mach numbers, with the angle the 
upstream magnetic field makes with the shock normal, Qsn , set to 5° . Solid lines are new results from this work using the Jiittner-Synge 
(J-S) equation of state and Fade' approximation described by Eq. ((26}. Dotted curves represent results from the previous work of Double et 
al. (2004). 



This approximation is accurate to 0.25% over the entire domain 
and is slightly less algebraically complicated than the approxi- 
mation in Eq. (38) of Double et al. (2004). By inverting Eq. (gSl) 
to obtain 'fgi , and using Eqs. and one can find 7gi 
function of r,; : 

- ^+R(r,)-1 ~ 3 + 10r, + 6rf ' ^^^> 

Inserting the approximation from Eq. (j26p in Eq. (|24l) provides 
Wi I Pi as a function of only , eliminating the fourth unknown 
quantity in Eqs. ([TT|) . ([T^ and p^ : this is the protocol adopted 
for the numerical solution of the Rankine-Hugoniot relations. 



3.4. Numerical Solutions for the Jump Conditions 

The resulting three equations were solved numerically using 
Mathematica. For the case of a plane-parallel shock, we com- 
pared directly with solutions displayed in Fig. lb of Heavens & 
Drury (1988) where downstream flow speeds are found as a func- 
tion of upstream flow speeds for parallel ( Qbh = 0° ) electron- 
positron shocks at various temperatures using the Jiittner-Synge 
equation of state just as we do. We find no observable differences 
between our results and theirs for plane-parallel shocks. For the 
case of oblique shocks, representative solutions, as functions of 
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Fig. 3. — Solutions to the Rankine-Hugoniot relations for r , GBt2 , and 0u2 , as in Figure (2] However, here, the upstream field obliquity is 
©Bfi ~ 85° , constituting quasi-perpendicular shocks. Values of the Alfvenic (A^a) and sonic (A^s) Mach numbers are as labelled. Again, 
solid lines are new results from this work using the Jiittner-Synge equation of state, while the dotted curves represent the solutions of Double 
et al. (2004). 



the upstream fluid rapidity Tij3i = Fis/^ixs : are displayed in 
Figs. [2] (quasi-parallel case) and |3] (quasi-perpendicular shock). 
The plots exhibit the velocity compression ratio r = /3ixs//32xs , 
and the downstream fluid frame field and fluid NIF velocity 
angles to the shock normal. Observe that hereafter, the sub- 
script "S" will be omitted when referring to NIF values for /3i 
and Fi . For the sake of comparison, these Figures also display 
equivalent plots for the same sonic and Alfvenic Mach numbers 
[as defined in Eq. (|4])] and shock obliquities, taken from Dou- 
ble et al. (2004). It is evident that the solutions here closely 
match those of Double et al. in both the non-relativistic and 
ultra-rclativistic regimes. The jump conditions reveal several 
characteristics, such as the declining compression ratio with de- 
clining Mach numbers of either variety, and small fluid defiec- 



tions at the shock in the ultra-relativistic regime. There are 
noticeable differences between our solutions and those of Dou- 
ble et al. (2004) in the trans-relativistic domain, but mainly for 
just the compression ratio. 

These differences in the two works, especially apparent for 
low sonic Mach numbers, are the result of two different assump- 
tions regarding the equation of state both upstream and down- 
stream of the shock. In this work, the effective jgi for a given 
flow speed and Mach number can be determined using equations 
(|?7|) and (g]). For A^s = 100,60, and 6, the upstream values 
of 7gi stay within 1% of the nominal non-relativistic value of 
5/3. However, in the Ais ~ 2.6 case, we find jgi « 1.6 for 
large values of Fi/Ji . In their work. Double et al. (2004) make 
an approximation assuming a cold upstream flow (i.e a large 
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sonic mach number with upstream 7^1 = 5/3 ), resulting in the 
following equation relating the downstream 7^2 to the down- 
stream flow speed rather than to the downstream pressure: 



7s2 



SlTrel - 1) 



2 

rel 



1 



i + 4r. 



3r 



rel 



where 



Prel 



1 - 



'1/2 



(28) 



(29) 



Note that the slight angle between the upstream and down- 
stream NIF flow velocity vectors spawns the approximation for 
r^e; ; the details and justification of this approximation can be 
found in section 3.1 of their paper. Assuming that 7^1 = 5/3 
then results in a small A% discrepancy in 7^2 relative to 
results from our Jiittner-Synge equation of state, in the lowest 
Mach number cases. From the plots, clearly the numerical eval- 
uations of the compression ratio are sensitive to the choice of 
the form of the downstream equation of state, i.e., 7^2 • 

In the limit that Fi/Si approaches infinity, the details of the 
equation of state become irrelevant, and a simple analytic solu- 
tion can be found to approximate the results of both approaches. 
Defining 



q = 2 

sin Qbh Pi 



Ml 
sin^ Ge 



(30) 



one can write the asymptotic limit for the compression ratio as 



Fi > 1 



(31) 



which corresponds to Eq. (47) of Double et al. (2004). The 
downstream fluid deflection and magnetic fields angles possess 
corresponding asymptotic forms for Fi 1 of 



tan 0Bf2 
tan 9^2s 



Fitanesfi \/r^~l , 

3-r 
2FitaneBii ' 



(32) 



When combined with Eq. ([9|), the asymptotic equation for 
tan6Bf2 becomes Eq. (40) of Double et al. (2004). Clearly, 
the fluid deflection is very small for ultra-relativistic flows, the 
hallmark of the extreme inertia of the upstream fluid. One can 
see that, for the range of sonic Mach numbers explored here, 
the critical parameter g is a comparatively weak function of 
the sonic Mach number inducing less than a 50% change in q 
as the sonic mach number varies from 2.6 to 60. In these cases, 
the Alfvenic Mach number and the upstream magnetic field an- 
gle &Bii are more important for determining the asymptotic 
behavior of the jump conditions. It is also evident that since 
r 3 when Fi ;3> 1 and g 3> 1 , the downstream fluid deflec- 
tion angle 6^23 in the shock frame is extremely small. 

4. RESULTS 

The simulation we have developed is capable of handling both 
subluminal and superluminal shocks of arbitrary obliquity. It 
can also simulate the effects of SAS or LAS with varying lev- 
els of cross-field diffusion controlled through the r] parameter. 
These broad capabilities encapsulate a range of properties that 
are relevant to astrophysical shock environs such as those in 



extragalactic jets in gamma-ray bursts and blazars. Moreover, 
they allow us to examine and expand upon a variety of previous 
investigations, including the semi-analytic studies of Kirk et al. 
(2000) and Kirk and Heavens (1989) as well as other simulations 
such as Ellison and Double (2004), and Niemiec and Ostrowski 
(2004). The following subsections highlight our key results in 
distinct parameter regimes: parallel shocks, oblique subluminal 
shocks, oblique superluminal shocks, and finally, LAS scenar- 
ios. Each of these sections provide physical scenarios where the 
power-law index is substantially different from the "canonical" 
a ~ —2.23 (where dn/dp = p~'^ ) result, which we demonstrate 
only holds at the shock location in parallel, ultra-relativistic 
(Fi ^ 1) shocks with a SAS scenario, concurring with previ- 
ous work. Altering this specific scenario yields power-law indices 
that depend on the microphysics (LAS vs. SAS), shock obliquity 
( ©Bfi ), and turbulence parameter (77) as well as the location 
relative to the shock front. A brief interpretation of these results 
in the context of blazars is offered in Section [5] 



4.1. Parallel Shocks 

Parallel shocks possess the important simplification that 
cross-field diffusion is immaterial. Accordingly, the model pa- 
rameter r] does not impact the spectra at the shock, and serves 
only to define the diffusive scale along the shock normal. For 
the case of relativistic parallel shocks, the canonical result a 
a ~ 2.23 power law spectrum was highlighted in the semi- 
analytic study of Kirk et al. (2000), but had been found previ- 
ously by Monte Carlo simulations (Bednarz & Ostrowski 1998; 
Baring 1999) and confirmed also by Ellison & Double (2004). 
However, exhibited results from these studies were spatially re- 
stricted to the immediate vicinity of the shock. The simulation 
presented here provides the opportunity to expand upon these 
studies and explore the spatial evolution of the particle distri- 
bution as well. The semi-analytical work of Kirk et al. (2000) 
provides the best basis for benchmarking simulated distributions 
at the shock discontinuity. Accordingly, in this study, shock pa- 
rameters are chosen in order to facilitate this comparison. The 
eigenfunction method of Kirk et al. (2000) built upon the ear- 
lier seminal work of Kirk & Schneider (1987) as an approach to 
solving the diffusion-convection equation in the neighborhood 
of a relativistic shock. Kirk et al. used this technique to gener- 
ate power-law indexes and angular distributions for accelerated 
particles at a strong, weakly-magnetized plane-parallel shock in 
the SAS limit. In this case, accelerated particles are defined 
as particles whose Lorentz factor far exceeds that of the shock, 
so that distribution characteristics in the injection domain are 
not traced. However, the injection process is modeled in Monte 
Carlo simulation approaches, and we show results probing this 
domain in later sections. 

Fig. m displays shock compression ratios, and our Monte Carlo 
results for spectral indices and angular distributions at the shock 
in the NIF, for parallel shocks spanning a range of rapidities 
Fi/3i . Moreover, it exhibits corresponding results from Kirk et 
al. (2000), and clearly illustrates that we find excellent agree- 
ment between Monte Carlo simulation results and Fig. 3 of their 
work. To facilitate comparison, we adopted the Jiittner-Synge 
EOS here, though we note that details of the shock parameters 
for Fig. 3 of Kirk et al. (2000) were not explicitly stated in their 
paper. This is, effectively, a case approximating that of the red 
curve in Fig. [2]here, save that Qbh — 0° . This minor change ac- 
tually simplifies the Rankine-Hugoniot solution and is shown as 
the solid black curve in left panel of Fig.|31 The spectral index a 
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Fig. 4. — Left Panel: The Rankinc-Hugoniot MHD compression ratio r and spectral indices a for noii-therinal particle distributions, as functions 
of the shock rapidity Vifji , for plane-parallel, 0Bfi = 0° , shocks. The compression ratio is computed according to the protocols of Section |3] The 
points for r and a correspond to select shock speeds with simulation data; most of these possess angular distributions illustrated in the right panel. The 
simulation runs were restricted to the small angle scattering (SAS) regime, for which the dashed curve corresponds to the low-magnetization semi-analytic 
determinations of a in Figure 4 of Kirk et al. (2000; labelled KGGAOO). Right Panel: NIF frame angular distribution functions, normalized to unity, 
for parallel (©Bfi = 0° ), relativistic shocks with bulk rapidities Fi/Si as labelled, and compression ratios r = 3.995,3.905,3.714,3.414 and 3.04, from 
lowest to highest speed (see points in left panel). Distributions were measured at the shock (x = 0) and sampled only high energy particles with rapidity 
7/9 ^ Fi/Si in each case. The simulation results are the histograms, acquired for the small angle scattering (SAS) case, and the smooth curves are the 
semi-analytic solutions that Kirk et al. (2000) obtained (see their Fig. 3) to the diffusion-convection equation in the SAS limit. 



is a monotonically-increasing function of TiPi , as in Kirk et al. 
(2000) and Baring (2004), reflecting the increased energization 
per shock crossing cycle that competes against the influence of 
a declining compression ratio. The angular distributions in the 
right panel of the figure closely match those from Fig. 3 of Kirk 
et al. (2000), all measured at the shock discontinuity. In this 
panel, < /i^ < 1 cases are for particles heading downstream, 
and < /is < are charges moving upstream. The distri- 
butions exhibit an increase in convective beaming downstream 
as the upstream flow speed increases. Such distributions were 
obtained as angle-dependent fluxes, and then divided through 
by the flux weighting factor /3fis before normalization. This 
introduces the apparent statistical noise in the neighborhood of 
/is = 0. 

To delve deeper into the anisotropics incurred in relativis- 
tic shocks, in Fig. [5] we examine the Fi/Ji = 10 case in more 
detail, extending the angular distribution illustrations to loca- 
tions upstream (left panel) and downstream (right panel) of the 
shock front. Again the distributions correspond to high energy 
particles with rapidity, 7/? » Fi/3i . In each panel, the black 
histogram is the distribution function at the shock, exhibited 
in Fig. m The origin of the shape of the angular distributions 
can be understood qualitatively. In non-relativistic shock sce- 
narios, high energy particles of speeds far in excess of that of the 
shock realize isotropy in all relevant frames of reference. How- 
ever, in relativistic shocks, even particles traveling very close 
to the speed of light can no longer be considered isotropic in 
all relevant frames and at all positions. Consider a relativistic 
particle returning to the upstream side of the shock from the 
downstream side. The upstream fluid frame velocity vector of 
this particle is initially pointed upstream. As the trajectory 



is perturbed via SAS seeded by field turbulence, the velocity 
vector in the upstream fluid frame performs a random walk. 
Because of relativistic beaming effects, once the particle's path 
is perturbed by the small angle 9 > 1/Fi in the upstream fluid 
frame, the shock frame x-component of the velocity (or angle 
cosine /i^ ) becomes positive, sweeping the particle back to the 
shock before it has the chance to isotropize in the upstream 
fluid frame. Accordingly, the parameter space around /i^ = 1 
is under-populated (actually exponentially suppressed) because 
the upstreaming particles have not had enough time to diffuse 
from fis < to fis > 0.9 before they are convected through the 
shock and downstream. This feature is critical to the hyper- 
efficient reflection in oblique relativistic shocks discussed in Sec. 
14.21 Note that at various non-zero obliquities, similar NIF frame 
angular distributions are elicited in the simulation at the shock, 
but are not shown. 

It is interesting to note that, as the particle detection plane 
is moved upstream, the domain of population suppression near 
/is = 1 expands to lower fig ■ This is because particles that have 
angle cosines closer to /x^ ~ — 1 are more likely to penetrate fur- 
ther upstream, so that when diffusing outside the Lorentz cone, 
they are less likely to populate near /i^ ~ 1 . This skews the dis- 
tribution progressively towards more negative /x^ . Additionally, 
the probability of particles reaching a position x upstream ex- 
ponentially declines with |x| on diffusive lengthscales (e.g. see 
Lee 1983, and Summerlin & Baring 2006, for illustrations of this 
in non-relativistic, heliospheric shock contexts). Accordingly, 
distribution functions taken at larger distances upstream suffer 
from poor statistics. Thus, the upstream distribution functions 
exhibited in the left panel of Fig. [5] are normalized to have the 
same area for display purposes. 
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Fig. 5. — Normalized NIF frame angular distribution functions for high energy particles with rapidity, 7/3 2> Fi/Ji , upstream (left panel) and down- 
stream (right panel) of the shock at various distances. The simulation run was for a parallel (©Bfi = 0° ), ultra-relativistic shock with Fi/Ji = 10 and 
compression ratio r = 3.04 , and diffusion in the SAS limit. Left panel: The black histogram is the distribution function at the shock and can be compared 
directly to the dotted line, which is an analytic result from Eq. (23) of Kirk ct al. (2000). The other 4 distribution functions are taken at increasingly 
large distances upstream of the shock. In units of Ti/Bimp/qB , the cyan curve is at x = —20 , the blue curve at a; = —80 , the green curve at x = —320 , 
and the red curve at x = —1280 . Distribution functions determined at larger distances upstream suffer from poor statistics, since few particles are able 
to diffuse so far upstream against the relativistic flow. Right panel: As in the left panel, the black histogram is the distribution function at the shock, 
with the other histograms now corresponding to a: = 400 (cyan), x = 1600 (blue), x = 6400 (green), and x = 25600 (red). The dashed line represents 
an isotropic distribution in the downstream fluid frame, as viewed by an observer in the NIF frame where the shock is stationary. 



The evolution of the distribution function downstream of the 
shock is shown in the right hand panel of Fig. [51 ranging from 
the distribution found at the shock (black histogram) to an 
isotropic distribution in the downstream plasma frame (red his- 
togram). As the particles move downstream, the relativistic 
beaming that biases the distribution to higher average values 
of is progressively enhanced, and they eventually isotropize 
in the downstream fluid frame (the red histogram in the right 
panel). The dashed line in that panel is the angular (density) 
distribution function in the shock frame, dNs/dfj,s , for particles 
that are isotropic in the fluid frame with a power-law distribu- 
tion dNf/dpf = 47rp^/(r, pf) oc pj"^ . Here /(r, pf) is the fluid 
frame phase space distribution function, which is a Lorentz in- 
variant. Hence in the shock rest frame, the angular distribution 

satisfies dNg/d^s oc plpf'^"'^'^^ for a fixed choice of Ps , which 
is imposed in this example by truncating the NIF distribution 
at a lower limit of pa = po . For ultrarclativistic particles, the 
relationship between ps and is given simply by the photon 
aberration formula pj = Ps^{^ — P^^s) ■ Accordingly, the angu- 
lar distribution in the NIF for isotropy in the fluid frame scales 
as dNg/dus c>c (1 — /?/.Js)^^'^+^) . This is what is illustrated in 
Fig. [5l which for the Fi ^ 1 case reduces to that for values 
/3 « 1/3 downstream and cr = 2.21 . 

This evolution of the anisotropy has consequences for the ob- 
served power-law index downstream of the shock. Because the 
average value of fj,s for the returning particles is lower than it 
would be for particles that are isotropizcd in the downstream 
fluid frame, the average bulk flow speed of the accelerated par- 
ticles is lower than U2x ■ As the angular distribution func- 



tion diffusively evolves toward isotropy in the downstream fluid 
frame at larger x , the average velocity of the particles also 
increases, asymptotically approaching the higher bulk velocity 
of the downstream thermal particles. This necessarily reduces 
the density of the high energy particles by conservation of par- 
ticle number flux. The scale length for the evolution of the 
distribution function is approximately the particle mean free 
path. Thus, higher energy particles with typically longer mean 
free paths isotropize farther downstream than do particles of 
lower energy. Accordingly, the spectral shape of the distribu- 
tion downstream of the shock evolves, illustrated in Fig. [6l in a 
manner that correlates with the spatial changes in the angular 
distribution. 

The black curve in Fig. [5] displays the distribution function 
at the shock. However, just downstream (cyan curve), the low 
energy particles have isotropizcd, thereby increasing their bulk 
flow speed and lowering their density, to conserve particle num- 
ber flux. Higher energy particles in this curve have not yet fully 
isotropizcd and possess slightly slower bulk speeds and thus, 
higher densities. So, if one is measuring distributions somewhat 
downstream of a relativistic shock, the observed power-law will 
be harder than the canonical a = —2.23 result realized exactly 
at X = . At large distances downstream of the shock, acceler- 
ated particles of all rapidities far in excess of V2P2X acquire the 
same beamed, isotropic distribution shown in the right panel 
of Fig. m so that their cumulative density adjustments during 
downstream diffusion are identical, and the power-law index re- 
turns to that realized at the shock. This variation of the power- 
law index with the position of a downstream observer relative to 
the relativistic shock represents a fundamental shift from non- 
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relativistic shocks, where the distribution function is isotropic in 
all relevant reference frames, and the spectral index is virtually 
independent of position when downstream of the shock. To our 
knowledge, this is the first time this effect in relativistic shocks 
has been highlighted in the literature. 




Fig. 6. — Accelerated particle distribution functions at positions up- 
stream and downstream of the shock. The positions and color coding 
corresponding to Fig. [5] with the histograms that fall sharply with de- 
creasing momentum corresponding to the upstream positions, and the his- 
tograms that make only a small adjustment from the distribution function 
at the shock corresponding to downstream positions. To clearly illustrate 
the differences between these distribution functions, the differential density 
distribution dN/dp has been multiplied by p^'^^ so as to generate zero 
power-law slope at the shock. The paucity of low energy particles at the 
upstream positions is due to their limited contraflow mobility, and is seen 
in non-relativistic shocks as the "convective peel-ofF' effect described in 
Summerlin & Baring (2006), see text for a discussion. The spectral vari- 
ations downstream of the shock are a result of energy-dependent density 
compression, and are addressed in the text. 



4.2. Oblique, Subluminal Shocks 

While small shifts in the power-law index can occur based 
on observation location and energy range in parallel shocks, the 
introduction of non-zero magnetic field obliquity creates more 
substantial ranges of power-law indices for the non-thermal par- 
ticle component. In non-relativistic shocks, the power-law index 
is independent of magnetic field obliquity, and depends only on 
the compression ratio: cr = (r + 2)/(r — 1) (Bell, 1978; Drury, 
1983; Jones & Ellison 1991). In relativistic shocks, the spectral 
index varies dramatically with obliquity, as will be exemplified 
in due course. In particular, the character of the spectral in- 
dex with respect to field obliquity hinges critically on whether 
the shock is subluminal or superluminal. Thus, our study of 
oblique relativistic shocks is divided into two sections to treat 
these parameter regimes separately. 

Consider first subluminal, oblique shocks in the small angle 
scattering limit. The first emphasis will be on the power-law 
behavior of the accelerated portion of the population; later on 
the injection efficiency will be addressed. In the simulation, 
for each run, the power-law regime is determined on an indi- 



vidual basis by inspection and can begin anywhere from 5 to 
100 times the mean injection (i.e. approximately downstream 
thermal) energy. A least squares fit in log-log space is used 
to determine the slope a . Results are depicted in Fig. [7] for 
/3i^ = ui^/c = 0.1,0.5 , and in Fig. E] for = 0.71 , for dif- 
ferent values of the turbulence or cross-field diffusion parameter 
77 = A/rg . The power- law index a is plotted as a function of the 
HT frame dimensionless speed /3iht = l^ix/ cosQbh ■ It is clear 
that there is a considerable range of indices a for non-thermal 
particles accelerated in mildly relativistic, oblique shocks. The 
essence of this array of indices and the global trends with Qbh 
and 77 were outlined briefly in Baring & Summerlin (2009) and 
Baring (2011), though a fuller interpretation ensues below. 

A feature of this plot is that the dependence of a on field 
obliquity is non-monotonic. When A/rg 3> 1 , the value of a 
at first declines as 0Bfi increases above zero, leading to very 
fiat spectra. As /3iht approaches and eventually exceeds unity, 
this trend reverses, and a then rapidly increases with increas- 
ing shock obliquity as the shocks become superluminal. This 
dramatic steepening is a consequence of inexorable convection 
of particles away downstream of the shock. The only way to 
ameliorate this rapid increase in index is to reduce 77 = X/vg to 
values below around 10 . Physically, this corresponds to increas- 
ing the hydromagnetic turbulence to high levels that force the 
particle diffusion to approach isotropy: = 1/(1 + 77^) in 

a kinetic theory description (e.g., Forman, Jokipii, and Owens 
(1974). This renders the field direction immaterial, and the 
shock behaves much like a parallel, subluminal shock in terms 
of its diffusive character. Note that this general character is also 
evinced in the very recent diffusion-convection equation analysis 
of Bell, Schure & Reville (2011) at shocks of lower speeds. Fig- 
ure 1 in their paper clearly illustrates that the distribution index 
hardens with increasing obliquity when the shock is well within 
the subluminal regime and softens when the luminal boundary 
cos &B11 = Pix is approached or crossed in quasi-perpendicular 
(and sometimes non-relativistic) shocks, unless the frequency of 
scattering is raised to the Bohm limit, for which the index then 
depends only weakly on the field obliquity. 

In studying this case, we again choose to use previous semi- 
analytic analyses as a benchmark for comparison: the work 
of Kirk and Heavens (1989, KII89 hereafter) is ideally suited 
for this purpose. KH89 calculated the spectral index of non- 
thermal particles at oblique, trans-relativistic shocks using the 
cigenfunction technique of Kirk & Schneider (1987) to solve the 
diffusion-convection equation. Their analysis was restricted to 
situations where particles do not diffuse across field lines; i.e., 
their collision operator contains only a pitch angle scattering 
term. They also assumed that particles conserve their mag- 
netic moment on crossing the shock, a standard analytic sim- 
plification. Results from Fig. 2 of their work are exhibited in 
Fig. [71 Note that their exploration was done exclusively in the 
HT frame and was thus limited to subluminal shocks. An inter- 
esting product of their work was the appearance of power-laws 
harder than both the non-relativistic and ultra-relativistic paral- 
lel shock results in Bcdnarz & Ostrowski, (1998) and Kirk et al. 
(2000), achieved as Uiht approaches the speed of light, but uix 
remains mildly-relativistic. This is an idealized result, because 
the limit of zero cross-field diffusion does not occur in Nature, 
since field turbulence abounds in astrophysical shocks, and is 
needed to drive acceleration. The Monte Carlo technique is ide- 
ally suited to examining how close to zero cross-field diffusion 
one must get to approach the particular analytical case explored 
by KH89. Fig. [7] clearly indicates that when X/vg ^ 10^ the 
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Fig. 7. — Power-law indices for simulation runs in the limit of small angle scattering (SAS), for an almost non-relativistic shock of upstream flow 
speed 0ix = uix/c = 0.1 {left panel), and a mildly-relativistic case with Pix = uix/c = 0.5 {right panel), for an MHD velocity compression ratio 
^ = Ulx/u2x = 4. The indices are displayed as functions of the effective de Hoffmann- Teller frame upstream flow speed /SijjT ~ /^ii/ cos 0Bfi , with 
select values of the fluid frame field obliquity Sefi marked at the top of the panel. Obliquities for which /3iht > 1 constitute superluminal shocks. The 
displayed simulation index results were obtained for different diffusive mean free paths A parallel to the mean field direction (see text), namely X/rg = 1 
(blue squares), A/r^ = 10 (black triangles), X/rg = 10^ (green pentagons), X/rg = 10'^ (red triangles) and X/rg = 10'* (magenta hexagons), as labelled. 



The lightweight black curve at the bottom labelled KH89 defines the semi-analytic result from Kirk < 
equation, corresponding to X/rg oo . 



: Heavens' (1989) solution to the diffusion-convection 



KH89 zero-cross field difFusion indices are closely reproduced for 
Pix = 0.1 and well approximated for /3ix = 0.5. The physical 
origin for these extremely hard power-laws will be discussed in 
Sec. HXTl 

To allow a direct comparison with the results of KH89, we 
adopted the same compression ratio of r = 4 , and the same for- 
mulation for the relationship between the upstream and down- 
stream magnetic fields. This formulation is found in equations 
(2), (3) and (4) of their work and is summarized below. It as- 
sumes a weak magnetic field that does not influence the plasma 
motion (i.e. Ma ^ 1 ); all the simulation runs used to generate 
Fig. [7] satisfied this high Alfvenic Mach number criterion. 



(33) 



Bo 



r?(r2-i)/3L 



1 



Given the upstream quantities above, and the compression ra- 
tio, r, we can solve for i?2HT and use our knowledge that B^ 
is constant across the shock along with an appropriate Lorentz 
transformation (see Sec. 13. 2p to find the appropriate down- 
stream value for B^ in any reference frame. 

In Fig. [71 while impressive agreement with the solutions of 
KH89 arises for /Si^^ = 0.1 when 77 ^ 10^ , for /3i^ = 0.5, 
we find that our simulation indices match the results of KH89 
at /3iHT <; 0.5 and just below /3iht = 1 , but are noticeably 
softer (higher a ) in the central part of the curve. Increasing 
rj as high as 10^ creates no appreciable change in the result- 
ing power-law index from that of 77 = lO'* : we believe we have 
reached the asymptotic limit of our simulation. Monte Carlo re- 



sults for (3ix = 0.3 are not depicted, but are similar to those for 
the ~ 0.5 case, and also exhibit a modest difference from 
the KII89 determinations of a at intermediate values of /3iht > 
while matching at the /3iht endpoints. We contend that the 
reason for the modest discrepancy between the two approaches 
probably lies in the assumption of conservation of magnetic mo- 
ment employed by KH89. This assumption facilitates an ana- 
lytic result but does not precisely describe orbiting particle re- 
flection and transmission properties at an oblique shock discon- 
tinuity (see Terasawa 1979; Drury 1983; Pesses & Decker 1986; 
for non-relativistic shock expositions). For parallel or quasi- 
perpendicular (in this case nearly-luminal) shocks, the mag- 
netic moment is conserved, and the two approaches converge. 
For obliquities in between, there is slight non-conservation of 
magnetic moment, and the precise tracing of gyro-orbits in the 
shock layer, as is enacted in the Monte Carlo technique, intro- 
duces modest but appreciable increases in a . 

It is imperative to go beyond the artificial r ~ 4 exploration, 
since relativistic shocks are somewhat weaker in their compres- 
sion. To this end we produced similar index plots for parame- 
ters more appropriate to internal shocks in GRBs and blazars 
using the Rankine-Hugoniot relations derived earlier. Specifi- 
cally, Fig. [5] displays Monte Carlo results for compression ratios 
that satisfy the Jiittner-Synge equation of state, Ais = 2.6 
( r = 3.02 ) and A^s = 60 ( r = 3.71 ), with the Alfvenic Mach 
number assumed large. The results mirror those in Fig. [7] in 
terms of overall character, with a large range of indices, cr ~ 1 
in near- luminal cases when X/rg > 10^ , and a rapid steepen- 
ing of the non-thermal distribution in superluminal cases unless 
X/vg <^ 10 . A particular index inferred from the radiation spec- 
trum of a single astronomical source can be accommodated by a 
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Fig. 8. — Power-law indices for simulation runs in the SAS limit, as in Fig. [T] but now for shock parameters more appropriate to the internal shocks 
associated with the relativistic jets that are believed to be the source of GRBs. Here, Pix = uix/c = 0.71 , and the compression ratio and sonic Mach num- 
ber are now r = U\xlu2x = 3.02 and A4s = 2.6 {left panel), and r = Uix/u2x = 3.71 and A^g = 60 {right panel), calculated via the Rankinc-Hugoniot 
relations derived in Sec. |3] The indices are again plotted versus the effective dc Hoffmann- Teller frame upstream flow speed Pinj^ = Pix/ cosOgfi , and 
selected fluid frame field obliquities ©Bfi are as marked at the top. The different diffusive mean free path cases were again A/rg = 1 (blue squares), 
A/rg = 10 (black triangles), A/rg = 10^ (green pentagons), A/rg = 10'^ (red triangles) and A/rg = lO** (magenta hexagons), as labelled. Also depicted 
are marker indices for three Fermi blazars, Mrk 421, 3C 66A and PKS 2155-304; they apply for arbitrary /3iht , and are truncated in the horizontal 
direction to aid clarity of the flgure. These mark the approximate expectation for a , uncertain to roughly ±0.2 , for an interpretation of the Fermi-hAT 
gamma-ray spectral indices as uncooled inverse Compton scattering {left panel) and strongly-cooled upscattering {right panel); see Sec.[5]for a discussion. 



range of choices for shock speed, Mach numbers, field obhquity, 
and turbulence parameter ry . This interpretative aspect is the 
subject of Section [SJ with a focus on blazars. 

Finally, note that Figs.[7]and[8]were prepared specifically with 
diffusion seeded by gyro-resonant interactions between charges 
and MHD turbulence in mind. In such cases, scattering descrip- 
tions are only physical if > 1 in Eq. ([T]), i.e. above the Bohm 
limit. Yet 77 < 1 regimes for diffusion can be realized for non- 
gyroresonant interactions with field turbulence that is perhaps 
grown via filamentation or Weibel instabilities. Trial simulation 
runs were performed in this 77 < 1 domain, and it was found 
that the distribution was not very sensitive to the choice of rj ; 
for example, reducing rj to 0.1 flattened the spectrum for the 
f3ix — 0.71 , /?iHT = 0.96 , r = 3.71 case by an index of around 
0.1 relative to that displayed in Fig.[S] This behavior is a conse- 
quence of diffusion in this "sub-Bohm" domain resembling that 
for the Bohm limit of 77 = 1 for gyroresonant diffusion. A more 
complete exploration of this domain is deferred to future work. 

4.2.1. The Action of Shock Drift Acceleration 

The bottom line of the preceding exposition is that the 
power-law index achieved in subluminal oblique shocks can 
be considerably less than even the a = {r + 2)/(r — 1) re- 
sult for non-relativistic shocks. For moderate obliquities and 
Tj = X/rg > 10^ it can become as hard as a = 1 . A power-law 
this hard can only be achieved in the case that particle escape 
from the acceleration region is miniscule. We illustrate here 
that the high rj cases that have low a are subject to strong 
shock drift acceleration (SDA), offering a close examination of 
the trajectories of energized particles that reveals prolonged re- 
tention in the acceleration process. A small fraction of particles 



incident from upstream can be reflected at the shock because 
they have suitable pitch angles, and these seed the retention in 
the acceleration process. A sample trajectory and associated 
momenta for a select particle undergoing such acceleration is 
displayed in Fig. [9] The particle was injected at super-thermal 
energies to circumvent improbable injection from the thermal 
population, a property that is discussed later in this subsection. 
The trajectory highlights two hallmarks of shock drift acceler- 
ation: coherent trapping in the shock layer, interspersed with 
upstream excursions after reflection at the shock (see Decker & 
Vlahos 1986, for illustrations in non-relativistic shock contexts). 
The reflection condition can be estimated by assuming conser- 
vation of magnetic moment in the HT frame, i.e., requiring that 
(1 - /^pi)/^iHT = (1 - Aip2)/-S2HT , where fipi is the particle's 
pitch angle cosine in the upstream (i=l) or downstream (i=2) 
HT frame. This assumption is technically valid only when par- 
ticles gyrate a large number of times in the shock layer (e.g., see 
Drury 1983), which arises when the gyroperiod is far inferior 
to the time it takes to convect one gyroradius downstream, i.e. 
when ps ^ Fi/JiHT ■ For magnetic moment conservation, given 
-B2HT > SiHT , it is clear that there are some values of fipi for 
which /ipi > cannot be satisfied. In these cases, particles are 
reflected rather than transmitted, and the shock is acting as a 
magnetic bottleneck. 

As particles gyrate in the shock layer, the work done dW on 
a charge can be computed using the Lorentz force, resulting in 
an equation 

{mcfl^ = p.^ = gp.|E+:^xB| ^ qp.E , (34) 
where E is the uxB drift field. In the uniform B fields either 
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Distance normal to shock, x/(r,/Sj^mc^/eBj) Momenta, (p,., Pg)/(r,|8,^mc) 

Fig. 9. — Left panel: A sample particle trajectory depicting strong SDA in an oblique mildly relativistic shock with dc Hoffmann- Teller frame speed, 
/3iHT = 0.96, NIF shock speed, /3ij; = 0.71, and compression ratio, r = 0ix/P2x = 3.71 in a low turbulence environment with X/vg = 10*. The 
projection is onto the x-y plane, where the u X B drifts lie in the y -direction. This portion illustrates two key features of shock drift acceleration: 
coherent trapping in the shock layer (colored red), interspersed with upstream excursions (black) after reflection via approximate conservation of magnetic 
moment. This particular particle gained orders of magnitude in energy before finding a pitch angle small enough to allow transmission and subsequent 
escape downstream. Right panel: The particle's position in the drift coordinate (y) direction, as a function of the magnitude of the momentum in the 
fluid frame, pf and the shock frame, ps ■ This illustrates a core property of shock drift acceleration: that over large times both pf and ps display 
approximately linear trends with the drift y . Prolonged energy gains occur only during shock layer gyrations, when the fluid frame momentum exhibits a 
"rectangular hysteresis" (red). The shock frame momentum possesses perturbed oscillatory temporal behavior (yellow) during the intervening upstream 
excursions. 



upstream or downstream, the energy gains and losses acquired 
during a gyroperiod exactly cancel, so that no net work is done, 
dW = mc^d'y = . In contrast, when a charge's gyromotion 
straddles the shock discontinuity, it samples the different elec- 
tromagnetic field on either side of the shock for different times, 
with the net acceleration on the upstream side of the shock 
being greater than the deceleration on the downstream side of 
the shock. Such an asymmetry in energy increments is explic- 
itly evident in Eq. (6) of Jokipii's (1982) exposition on SDA 
in non-relativistic shocks (see also Webb, Axford & Terasawa 
1983, for relativistic cases). This energization can be seen in 
the counter-clockwise rotations of the p/ curve in Fig. [HI In 
this curve, vertical motion indicates travel upstream or down- 
stream of the shock where Pf is constant. The horizontal lines 
are shock crossings where the local pf changes in transits of the 
shock discontinuity. By relating elapsed time during SDA to in- 
crements in the drift coordinate y , it is simply shown (Jokipii 
1982) that 

dW = mc^dj = qEydy , (35) 

or d{p/mc) oc [eBi/mc^) dy in the relativistic limit of 7 ^ 1 . 
This energy gain scales linearly with displacement along y , with 
the proportionality constant depending on the shock obliquity. 
Such a linear scaling (in ps ) is the punchline of the right panel 
of Fig. [HI where y effectively represents a time coordinate dur- 
ing shock drift episodes. This energization proportionality is 
clearly evident during shock interactions for the selected par- 
ticle, and is an established hallmark of SDA. But notably, on 
larger scales, for both pf and Ps , it is also an approximate de- 
scription of the cumulative energization spanning multiple shock 
encounters. This follows as a consequence of the constancy of 



Pf and only moderate changes to y during upstream excur- 
sions; i.e. material energy changes arise only during shock in- 
teractions. The Ps curve's behavior shows smooth and rapid 
energization during gyrations in the shock layer followed by a 
quasi-oscillatory epoch for ps coupled with random spatial dif- 
fusion during the upstream excursions. The gradual increase 
of Ps during the upstream excursions results from the angular 
diffusion of the particle's fluid frame momentum vector in the 
upstream region. This increases Ps as the fluid frame velocity 
vector becomes more aligned with the flow vector. The widening 
of the gyrations in y is also due to angular diffusion gradually 
increasing the pitch angle upstream of the shock. 

In typical SDA, a reflected particle gains energy and is sent 
back upstream to encounter the shock again. However, this 
alone can not create such an unusually hard power-law. In 
the case of non-relativistic shocks, this process happens com- 
monly in oblique shocks when low levels of turbulence are as- 
sumed, with no impact on the power-law index: it still retains 
its cr = (r + 2)/(r - 1) value (sec Jokipii 1982; Lee 1984; Arm- 
strong et al., 1985; Decker and Vlahos 1986; Webb, Axford & 
Teresawa 1983; Pesses et al, 1982; Decker, 1988; Vandas, 2001). 
This is because the energy gained through shock drift acceler- 
ation (SDA) is exactly canceled by decreases in the efficiency 
of diffusive shock acceleration (DSA) in oblique shocks, since 
the angular distribution is effectively isotropic in the NIF and 
fluid frames. Furthermore, this isotropy guarantees that select 
particles arc not trapped in the shock layer for long durations, 
because they have the same probability of being reflected at 
each shock encounter. However, in the case of oblique sublumi- 
nal shocks whose HT frame speed is approaching the speed of 
light, we find different circumstances. 
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As discussed in section 14.11 in relativistic, parallel shocks in 
an SAS scenario, the field-aligned component of the distribu- 
tion function is depicted at the shock. While the first shock en- 
counter is that of an isotropic plasma boosted by the upstream 
flow speed with a strongly field-aligned population of particles, 
the distribution of particles returning after being reflected the 
first time is missing these field aligned particles because the 
population has, in general, had insufficient time for any of them 
to become field aligned in the downstream direction before re- 
encountering the shock. This is a result of the fact that rela- 
tivistic shocks preclude particle speeds from far exceeding the 
relevant flow speeds and becoming effectively isotropic in all ref- 
erence frames. Subluminal, oblique relativistic shocks exhibit 
this same behavior. Upstreaming particles must initially have 
pitch angle /Xht ~ — 1 in order to travel upstream faster than 
the field line they are on is being convected downstream. Once 
they diffuse in pitch angle, they are quickly swept back into the 
shock before they achieve isotropy in the HT frame. Thus, as 
in the case of parallel shocks, the field aligned, /Iht ~ 1 , com- 
ponent of the angular distribution function is depleted at the 
shock. Yet only such field-aligned particles are capable of pene- 
trating the magnetic bottleneck at the shock. The net result of 
this dearth of field-aligned particles is a significant reduction in 
the fraction of particles transmitted through the shock at each 
encounter. This enhances the probability of reflection harden- 
ing the power-law. As the HT frame speed approaches c, the 
transmission region can become completely depopulated leading 
to near 100% reflection and the observed a = 1 power-law. 

The value of this index value is dictated by energy arguments. 
For long-lived trapping of select particles in the shock layer, en- 
ergy in the particle population is transmitted from one Lorentz 
factor bin [7, 7-1-^7] to the next one above, with minisculc loss. 
The energy content of this bin is jN , and when deposited in the 
next bin above, it is increased by SDA by an amount oc dy (x dj 
according to Eq. (|35|) , which is independent of the value of 7 , 
but just on the electromagnetic structure of the shock layer. It 
then follows that the energy increment in going from adjacent 
bins is 7c?A^ cx ^7 , so that dN/ dj (x I/7 , i.e. cr = 1 . Introduc- 
ing significant losses reduces this energy increment, and thereby 
steepens the spectrum. Along with the offering in Baring & 
Summerlin (2009), this is the first identification of the impor- 
tant role shock drift acceleration can play in determining the 
spectral index in relativistic shocks. In non-relativistic shocks, 
SDA does not influence the spectral index. 

Note that the highly enhanced action of SDA is restricted to 
high 77 and SAS regimes. Invoking a LAS scenario completely 
eliminates the enhanced probability of reflection as particles are 
isotropic after their first upstream scattering and encounter the 
shock as such: all subsequent shock encounters have the same 
probability of reflection as the flrst shock encounter. For SAS 
scattering, increasing the amount of turbulence present (reduc- 
ing r] ) allows particles to scatter into the transmission cone and 
subsequently softens the power-law as shown in Figs. [7] and [8] 
Trajectories for such rj < 10^ cases (not shown) exhibit a more 
"wonky" gyration and reveal prompt convection downstream, 
shutting down the opportunity for repeated episodes of coher- 
ent acceleration. 

A few comments are necessary on the feasibility of encounter- 
ing parameters that drive hyper-efficient SDA. The extremely 
large values of rj required correspond to very low levels of tur- 
bulence that are not anticipated near shocks: it would require 
a truly unusual set of physical parameters to produce power 
laws significantly harder than cr = 1.3 using this mechanism. It 



should also be noted that the bulk thermal particles that create 
these strong shocks must be cold compared to the flow speed in 
order for the shock to form. Thus, despite the fact that they 
receive substantial kinetic heating during their flrst shock cross- 
ing, they will not meet the p ^ FiUiht criterion and will have a 
reduced chance for reflection. Particles may have only a couple 
gyro-orbits that pass through the shock during an encounter and 
the range of phases that permit them to reflect is reduced dra- 
matically. For V < UiHT , it is physically impossible for particles 
to diffuse upstream along fleld lines and some amount of cross- 
fleld diffusion is necessary for particles to return to the shock 
at all. This creates an injection problem; see Ellison, Baring & 
Jones (1995) for a discussion of this in non-relativistic, oblique 
shock contexts. 

One naturally asks how these energetic particles that so effi- 
ciently participate in SDA get accelerated to high energies in the 
flrst place. Fig. [TUl illustrates the injection problem for r = 3.71 
shocks with very warm particles ( A^s = 4). The distributions 
in these plots were used to determine the spectral slopes dis- 
played in Fig. [S] for each value of /3iht ■ High turbulence envi- 
ronments are clearly able to inject the particles efficiently, but 
the power-law is steeper, namely a is higher. Low turbulence 
environments have almost no injection until particles achieve 
V > MiHT , at which point a strong, low a power-law devel- 
ops as particles become trapped in the shock drift mechanism. 
This becomes particularly pronounced in the /3iht = 0.96 case, 
an almost luminal shock situation, where convection of thermal 
upstream particles through and downstream of the shock is ex- 
tremely rapid. The rapid decline of injection efficiency with rj 
is an important factor in interpreting the action of shock accel- 
eration in astronomical sources, an issue discussed in Section [S] 

4.3. Oblique, Superluminal Shocks 

In superluminal shocks, it is physically impossible for even 
particles moving at the speed of light to diffuse upstream along 
the field lines. Without any cross-field diffusion, the particles are 
inexorably swept downstream after passing through the shock 
without reflection, regardless of their energy; see Begelman & 
Kirk (1990) for illustration of such trajectories. Thus, the parti- 
cles in the low turbulence environment that relied on reflection 
for retention in the system are lost and not accelerated. Con- 
sequently, the power-law becomes very soft. Particles in high 
turbulence environments are not truly affected by the change 
in the obliquity of the magnetic fleld because they can travel 
across fleld lines as easily as they can along them. It is gen- 
erally found that for /3iht < 1 decreased turbulence enhances 
acceleration in the power-law tail, and for /3iht > 1 , the op- 
posite holds. This result adds new perspective to the paper 
by Ellison and Double (2004), which presented results showing 
that the power-law tail indices in ultra-relativistic ( Fi ^ 1 ) 
shocks are extremely sensitive to both t] and the obliquity of 
the magnetic field, with the power-law index increasing sharply 
as these parameters increase. These dependences are also seen 
in the mildly-relativistic shocks discussed here, but with some- 
what less sensitivity to rj and Qbh ■ 

While it is clear that an increase in rj will soften the power- 
law in oblique, superluminal shocks, the exact values of the 
power-law index are simulation-dependent. Therefore, it is pru- 
dent to compare our results with those of Ellison and Double 
(2004, hereafter ED04). The code used in ED04 is a Monte Carlo 
simulation that is algorithmically very similar to the simulation 
presented in this paper. However, the simulations were devel- 
oped independently and can each serve as an objective test of 



19 



0.71, /SjH^ = 0.84 




1 2 
Logio[p/(mc)] 



T3 



o 



O 



o 




-14 



1 2 
Logio[p/{mc)] 



Fig. 10. — Full particle distributions for simulation runs in the small angle scattering (SAS) limit, for the strong mildly-relativistic shocks of upstream 
flow speed f^ix = uix/c = 0.71 whose indices are displayed in Fig. [8] Here the de Hoffmann- Teller frame upstream flow speed was set at /^ifjT ~ 0.75 
( ©Bfi ~ 48.2°), and the five values of the diffusive mean fee path A/r^ = 1, 10, 10^, 10'^, 10^ correspond to those in Fig.|8]- the color-coding of the 
distributions and the spectral index results in the respective Figures coincides. The velocity compression ratio was fixed at r = Uix/u2x = 3.71 , and the 
upstream temperature corresponded to a sonic Mach number of Mg = 4.04 ( T = 5.45 X 10^ K for an e~ - e+ shock). 



the other. For non-relativistic shocks, both simulations find the 
standard results of cr = (r + 2)/(r — 1) where r is the compres- 
sion ratio of the shock and a is defined such that dn/dp ~ p^" ■ 
In the case of ultra-relativistic, parallel shocks, both simula- 
tions also find the theoretical results a = —2.23. In the regime 
of oblique, relativistic shocks, ED04 focused predominantly on 
superluminal cases of high Mach numbers, with high levels of 
turbulence near the Bohm diffusion limit. 

In Fig. [TT] we compare results from our simulation (his- 
tograms) to the results from the top panel of Fig. 7 in ED04 
(solid lines) for a relativistic shock with upstream Fi = 10 
and compression ratio, r = 3.02 , for different values of the up- 
stream magnetic field obliquity, Gen , and 77 = A/r^ . Both 
sets of results are in the SAS limit, with the Rankine-Hugoniot 
solutions for the MHD shock obtained using the prescription 
of Double et al. (2004), as outlined in ED04, as opposed to 
the Jiittner-Synge EOS scenario; see Sec. 13.31 for details. For 
©Bfi = 0° , both simulations produce a result very close to the 
canonical cr = —2.23 power-law (note that the y-axis is multi- 
plied by p2.23 gycj^ a horizontal line is the canonical result). 
However, at 9Bti = 60°, ui^^l cos Gen ~ 2 making the shock de- 
cidedly superluminal. Thus, cross-field diffusion is essential in 
order for particles to be able to return to the shock and increas- 
ing T] increases the power-law index considerably. This trend 
is identified by both simulations in this ?y < 6 domain; in our 
simulation it continues to somewhat higher t] before statistical 
degradation inhibits determination of the spectrum. In general 
character, results from the two simulations are clearly similar, 
and numerically they are close. Yet there is a real difference, 
with the index determination differing between the two sets of 
results within the range of 1 — 2%. The numerical precision 
of (T in our simulation is of the order of 1%. Comparisons 
were made with other superluminal shock results published in 
ED04, finding similar levels of agreement. The origin of the 



small spectral differences evinced in Fig.[TT]is not yet clear; we 
believe that the thorough comparison of our indices and angu- 
lar distributions with the semi-analytic approaches of Kirk & 
Heavens (1989) and Kirk et al. (2000), among other simulation 
checks, advocates for the robustness of our results. 

The simulation developed by Niemiec and Ostrowski (2004, 
hereafter NO04; see also Niemiec & Ostrowski 2006) provides 
another opportunity for comparison. NO04 uses a fundamen- 
tally different mechanism for particle scattering from our Monte 
Carlo simulation. Instead of phenomenologically scattering par- 
ticles and specifying a momentum dependence for the mean free 
path, their simulation injects a prescribed spectrum of turbulent 
magnetic field structure that is superposed on the bulk magnetic 
field. Variations in the magnetic field perturb the gyro-orbits of 
the particles, and, are intended to mimic turbulence that a par- 
ticle might encounter. Because NO04 uses a turbulent magnetic 
field, particle trajectories must be integrated over much shorter 
time steps than is possible in the Monte Carlo code presented 
in this paper. This necessarily results in longer run times and 
poorer statistics. These poor statistics are particularly evident 
in the angular distributions produced in Figs. 4, 8, and 11 of 
their work. We find qualitatively similar general behavior for 
results from the two techniques. 

Consider Fig. 2 of NO04, where distribution functions for 
accelerated particles in subluminal, mildly-relativistic oblique 
shocks with compression ratio, r — 5.11 are exhibited. This 
value of r exceeds the non-relativistic, strong shock limit of 
r = 4 , and is appropriate for the mildly relativistic electron- 
proton shocks studied in Heavens & Drury (1988), from whom 
they acquired their compression ratios. Though there is not 
a simple relationship between the two turbulence parameters, 
{5B/B in their work, and r] in our simulation) they are cor- 
related with low values of 5B/B corresponding to high values 
of T] . In subluminal shocks for small 5B / B , the power law 
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index for the distribution function, f{p), which they call a, is 
approximately 3. This is related to the power-law index for the 
differential density distribution by cr = a — 2 or cr « 1 in their 
low turbulence limit. This agrees with both our results and 
those of KH89. Additionally large amounts of turbulence soften 
the power-law as we observed in our simulations runs. Finally, 
it is interesting to note that at high energies, where particles are 
resonant with larger wavelengths than the stirring scale of their 
simulation, turbulence disappears for particles at these energies 
(?7 — ^ oo) and the power-law becomes a ~ 1 , in good agree- 
ment with our results. In superluminal oblique shocks such as 
those shown in Fig. 5 of their paper, although the statistics are 
somewhat poor, it is clear that low turbulence is no longer an 
asset to the acceleration process. In the shocks shown there, the 
power-law is consistently softer in the low turbulence cases, and 
in some cases no acceleration occurs without significant turbu- 
lence. In the high turbulence case, power-laws are still produced, 
but then the definition of the shock obliquity and the details of 
particle diffusion become more important in determining the re- 
sulting power-law index, rendering comparison with our results 
less insightful. 
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Fig. 11. — A direct comparison between distribution results from our 
simulation, the histograms, and those Ellison and Double (2004), the solid 
lines, for a relativistic shock with upstream Lorcntz factor Fi = 10 and 
compression ratio, r = 3.02, for different values of the upstream magnetic 
field obliquity, 0Bfi , and rj = X/vg , the turbulence parameter. Specifi- 
cally, the results from the top panel of their Fig. 7 are compared here with 
the y-axis multiplied by p4.23 |-q snatch the presentation in that figure. 
Both simulations identify the same trends with only minor differences in 
the value of the power-law slope. 

4.4. Large Angle Scattering Domains 

One aspect of the simulation parameter space that has been 
neglected until now is the impact of varying the microphysics of 
the turbulent interactions; all previous results have focused on 
the small angle scattering limit. In this section, we explore such 
using our Monte Carlo simulation to model relativistic parallel 
shocks, by varying 6'scatt , the angular width of the conical sec- 
tor into which the particle's momentum vector is scattered at 
each encounter with magnetic turbulence. A value of ^scatt = tt 



corresponds to large angle scattering (LAS), where the particles 
scatter ~ 1 time per mean free path; this is the domain first 
highlighted by Ellison et al. (1990a). A small value corresponds 
to SAS, where the particles scatter N times per mean free path, 
where N is given by Eq. ©. 

Fig. [12] depicts accelerated particle distributions for two dif- 
ferent shock speeds, illustrating the multitude of power-law in- 
dices available while varying only the scattering angle, ^scatt , 
and fixing the obliquity at 9Bfi = 0° . Observe that for such 
parallel shocks, the distributions are independent of the dif- 
fusion parameter 77 since they are measured just downstream 
of the shock. This depiction complements results published in 
Fig. 2 of Stecker, Baring & Summerlin (2007), and illustrates 
two primary results. The first is that LAS scattering produces 
a step-like structure in the accelerated distribution, a character- 
istic first identified by Ellison et al. (1990a). Each step corre- 
sponding to particles with increasing numbers of shock transits. 
In other words, the first step consists almost entirely of particles 
that have crossed the shock 3 times. Particles in the second step 
have almost all crossed the shock 5 times, etc.; see Baring (2004) 
for an illustration of this correlation. The precise correlation be- 
tween shock transit number and particle energy weakens as the 
structure damps into a power-law. The prevalence of the step 
structure, and how high in energy it extends before relaxing into 
a power-law, increases with the Lorentz factor Fi of the shock. 
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Fig. 12. — Particle distribution functions dN/dp from parallel shocks 
that are either mildly-relativistic {Fi/3i = 3, i.e. 01 = ui/c 0.949) or 
ultra-rclativistic ( Fi/3i = 10 , i.e. /3i = ui/c 0.995 ; multiplied by 10^ to 
effect clarity of depiction), of velocity compression ratios r = ui/u2 ~ 3.24 
and r ~ 3 , respectively. For these simulation runs, scattering off hydro- 
magnetic turbulence was modeled by randomly dcficcting particle momenta 
by an angle within a cone of half-angle Sgcatt 1 whose axis coincides with 
the particle momentum prior to scattering. Values of ^?scatt span the range 
from large angle scattering (LAS: fgcatt < t 3> 1/Fi ) to small angle scat- 
tering or pitch angle diffusion, when fgcatt < l/^i and the distributions 
become independent of the choice of Sg^att • AH distributions asymptoti- 
cally approach power-laws dN/dp oc p~'^ at high energies. For two LAS 
cases, these power-laws are indicated by lightweight lines, with indices of 
(T = 1.61 (Fi/?! = 3) and cr = 1.62 (Fi/3i = 10). 
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The second major result is that decreasing the scattering an- 
gle removes this structure but at the same time, softens the re- 
sulting power-law. A complete investigation of why the power- 
law is harder in the LAS scenario is deferred to future work, 
yet the origin of this trend in ^scatt centers on the distribution 
function at the shock. While LAS produces a beamed isotropic 
distribution similar to the red histogram of the right panel of 
Fig. [S]at the shock, SAS generates a distribution like that shown 
in the black histogram of the same panel. Both the probabil- 
ity of return to the shock from the downstream side, and the 
energization per shock crossing, are functions of this angular 
distribution function (Bell, 1978; Peacock, 1981; Blasi & Vi- 
etri 2005). For non-relativistic shocks, the distribution func- 
tion is necessarily isotropic to leading order, which restricts the 
power-law index to be only a function of the compression ra- 
tio. Relativistic beaming is the probable cause for breaking this 
degeneracy in ^scatt space. Evolution of the angular distribu- 
tion with ^scatt can be inferred from Fig. 6 of Blasi & Vietri 
(2005). Note that introducing magnetic field obliquity can alter 
the nature of this trend, as is indicated in Morlino, Blasi & Vi- 
etri (2009). Future work will explore the variation of asymptotic 
values of the index a as functions of ^scatt and Osfi , and also 
compare the Monte Carlo simulation values with those obtained 
from the semi-analytic, transport equation approach of Blasi & 
Vietri (2005) and Morlino, Blasi & Vietri (2009). 

5. OBSERVATIONAL CONNECTION: THE IMPACT ON BLAZAR 
GAMMA-RAY INTERPRETATION 

To briefly outline how these simulation results are relevant to 
astrophysical contexts, we discuss blazars, the subset of active 
galactic nuclei possessing relativistic jets of material emanating 
from the supermassive black holes at their centers; these jets 
are oriented virtually towards the observer. Blazars were dis- 
covered as a class of gamma-ray sources by the EGRET experi- 
ment on the Compton Gamma-Ray Observatory (Hartmann et 
al. 1992), and subsequently detected by ground-based Qercnkov 
telescopes at TeV energies (Punch et al. 1992). The EGRET 
blazar measurements have been built upon in the last three years 
by Fermi Large-Area Telescope (LAT) detections of dozens of 
blazars, offering improved spectroscopy. The TeV-band signals 
typically exhibit steep photon spectra (e.g. see Krennrich et al. 
2002; and Aharonian et al. 2003, for observations of Mrk 421) 
that include the absorption due to pair producing interactions 
77 e^e~ with infra-red and optical light generated by the 
intergalactic medium along the line of sight to the observer. Ex- 
tremely flat particle distributions are inferred in some blazars 
after correcting for this attenuation (see, for example, Stecker, 
Baring & Summerlin 2007), with indices as low as cr ^ 1.5 in 
high rcdshift sources. Coupled with the TeV-band capability, 
the Jermz-LAT detections of blazars enable refined diagnostics 
by extending the observational window over a much larger en- 
ergy range, and most crucially, including below the 77 — > e^e~ 
attenuation window. Accordingly, Fermi observations can probe 
more directly the underlying radiating particle population. The 
implications of this we explore here. The reader can consult 
Baring (2011) and references therein for the interpretation of 
relativistic shock acceleration in gamma-ray burst contexts. 

Pertinent blazar data from the Fermi-LAT and TcV telescopes 
can be found in the GeV-TeV blazar "compendium" in Abdo 
et al. (2009). There is also the more extensive AGN catalog 
of Fermi in Abdo et al. (2010). For the purposes of discern- 
ing indices a of particle populations generating the gamma-ray 
emission, it is important to consider photon spectra below any 



turnovers that may appear in the LAT band. This biases the 
data selection to below 1 GeV, and a nice tabulation of this 
for Fermi-LAT blazars is given in Tables 5 and 6 of Abdo et 
al. (2009). Therein, and in the various spectral plots given in 
that paper, it is clear that there is a modest spectral steepening 
above around 1 GeV in around 50% of Fermi-LAT blazars; this 
becomes much more pronounced above 100 GeV. From this data 
compilation, we use photon indices of = 1.72 for PKS 2155- 
304, = 1.78 for Mrk 421, and = 1.97 for 3C 66A, as a 
sample blazar selection. The uncertainties on these indices are 
of the order of ±0.1 , which propagate into inferred particle in- 
dices a . Note that a sizable fraction of LAT band AGN indices 
in the catalog of Abdo et al. (2010) fall below 2 , so that 

our choice here is reasonable. Note also that 3C 279 is consid- 
erably steeper in the LAT window, which could be a signature 
of a low-energy onset of the spectral turnover, or the operation 
of the Klein-Nishina regime of inverse Compton scattering. 

Consider first a standard Icptonic model interpretation of 
blazar emission as emanating from inverse Compton scattering 
by shock-accelerated electrons upscattering low energy photons. 
The relationship between the particle index a and the photon 
one is then a = 2a^ — 1 (e.g. see Chapter 7 of Rybicki & Light- 
man 1979), if there is insignificant radiational cooling. This 
applies to synchrotron self-Compton (SSC) scenarios, where a 
single population of electrons emits the synchrotron radiation 
that it then upscatters to the gamma-ray band, or to an exter- 
nal supply of seed photons. For our select blazars, these indices 
fall in the range 2.44 < a < 2.94 , and are marked on the left 
panel of Fig. |8l They have an uncertainty Act ~ ±0.2 . From 
this it is clear that for shock-layer diffusive scattering in the 
SAS regime, acceleration at mildly superluminal oblique shocks 
provides a good description for all three of the blazars. This re- 
quires strong scattering, X/rg < 10 . If X/rg ^ 2 , i.e., near the 
Bohm diffusion limit, then Mrk"421 and PKS 2155-304 could be 
modeled with subluminal shocks. 

To contrast this, consider an alternative picture of "strong 
cooling" by inverse Compton scattering (or synchrotron radia- 
tion). This corresponds to rapid acceleration of the leptons at 
shocks, followed by convection and diffusion away downstream 
into a larger radiative zone where the gamma-ray signal is gen- 
erated over longer timescales. Then, as is well known, the time- 
averaged effective electron distribution that radiates is a power- 
law of index a+l (e.g., see Blumenthal, 1971, for an analysis). 
The steepening reflects a pile-up at lower Lorentz factors 7e in- 
duced by the fact that the energy loss rate for both synchrotron 
and inverse Compton cooling of electrons scales as 7g . The 
consequence is that now the relationship between the particle 
index a and the photon one is ct = 2a-y ~ 2 . These indices 
are marked on the right panel of Fig. |S1 and exhibit a shift 
of unity down from those in the uncooled case. Again, they 
have an uncertainty Act ^ ±0.2 . Now, subluminal regimes 
are clearly suggested, and in the cases of Mrk 421 and PKS 
2155-304, move the inferences into 77 = X/rg > 10 territory. 
This opens up the question of whether such weak turbulence 
can persist in a blazar jet shock. However, there is considerable 
uncertainty in the spectral data, and increasing the scattering 
angle ^scatt somewhat above the SAS limit will reduce the value 
of 77 required to generate a particular value of a . 

These inferences should be viewed only as general guidelines, 
modulo the uncertainties in ct spawned by the precision of 
Fermi-LAT spectral index determination. We have selected one 
particular shock speed and restricted the discussion to SAS 
regimes. Clearly, there is a range of shock speeds, field obliq- 
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uities, scattering angles ^scatt and turbulence parameters rj 
that can satisfy a measured gamma-ray index. In addition, this 
discussion has focused on leptonic models; in hadronic ones me- 
diated by pion decay, generally a-y approximately traces the 
particle index a , mimicking the situation in the right panel 
of Fig. m Significant Klein-Nishina modifications to the spec- 
trum in the LAT window can further complicate the interpre- 
tation. While more extensive study can hone the parameter 
space, broader gamma-ray coverage below 100 MeV and multi- 
wavelength modeling are necessary to make significant strides. 
The multiwavelength aspect is more immediate in terms of its 
possibilities. For example, if the acceleration at the shock is lim- 
ited by synchrotron cooling in controlling the maximum electron 
Lorentz factor, then the turnover in the synchrotron component 
is at an energy of ~ nieC^ / {afr]) (for af = e?/hc as the fine 
structure constant), a well-known result that is discussed in Gar- 
son, Baring & Krawczynski (2010) with reference to Mrk 421. 
Evidently, this energy must match that seen in hard X-ray/soft 
gamma-ray observations, providing additional constraints on ?/ 
that often may not be near the Bohm limit of 77 = 1 . This il- 
lustration serves to motivate future multiwavelength models of 
blazar spectra using complete distribution functions from accel- 
eration simulations like those presented in this paper. 

6. CONCLUSIONS 

This paper has presented new results from a robust Monte 
Carlo simulation that complement and extend previous semi- 
analytic and computational results. It employs the simulation 
technique devised by Ellison et al. (1981) that was extended to 
relativistic shocks by Ellison et al. (1990a). The simulation pro- 
duces steady state distribution functions for planar shocks of in- 
finite extent for large ranges of shock speeds, energies, and posi- 
tions, simulating both the injection and acceleration of particles 
via first-order diffusive shock acceleration. By using the unique 
advantages that a simulation has over semi-analytic, diffusion- 
convection equation solution techniques, we are able to expand 
upon the work of previous authors by examining various turbu- 
lence regimes and probing individual particle trajectories. This 
affords specific insights that cannot be gleaned from idealized 
cases that are analytically tractable. Our body of results leads 
to several key conclusions: 

• The power-law index in relativistic shocks samples a con- 
siderable range of values, and depends critically on the na- 
ture and magnitude of turbulence, the shock speed, and 
the shock field obliquity. This range extends from ex- 
tremely hard power-laws with cr « 1 to extremely steep 
distributions where simulation statistics preclude discern- 
ment of acceleration beyond the thermal injection domain. 
Notably, ultra-relativistic shocks do not necessarily pos- 
sess the "canonical" a k, 2.23 power-law distribution, a 
result evident in the previous works of Kirk & Heavens 
(1989), Ellison & Double (2004), and Stecker, Baring & 
Summerlin (2007). 

• When small angle scattering (SAS) is invoked, the value 
of uix/ cos 0Bfi defines a critical division point in the pa- 
rameter space. When it is less than c , oblique shocks 
in low levels of turbulence accelerate high energy parti- 
cles extremely efficiently via shock drift acceleration, but 
when this quantity is greater than c , turbulence becomes 
vital to injection and acceleration in oblique, superluminal 
shocks. In both cases, weak levels of turbulence strongly 
inhibit injection from the thermal population. 



• Invoking large angle scattering produces significant struc- 
ture in the high energy particle distributions in relativistic 
shocks, a phenomenon first identified by Ellison, Jones & 
Reynolds (1990), but also generates slightly harder distri- 
butions than a similar shock in the SAS scenario, where 
there is only a power-law with little discernible structure. 

These results represent important advances for determining the 
nature of turbulent shock environs in blazar and gamma-ray 
burst jets, and in other astrophysical objects. Such interpre- 
tations are, admittedly, complicated by the particular spatial 
environment and radiation emission mechanism chosen for gen- 
erating the observed photon spectra from these sources. Yet 
global insights such as deciding between subluminal or superlu- 
minal shock environments are now possible. 

To develop our model to aid future interpretations of astro- 
physical shocks, additional details of shock physics will be incor- 
porated into the simulation. In shocks such as those discussed 
above with cr < 2 , a majority of the energy in the system will be 
found in the accelerated particles. The Rankine-Hugoniot jump 
conditions for the shock must then be modified, since a step 
function shock profile is no longer a valid approximation. This 
yields a non-linear acceleration phenomenon that is already seen 
clearly in the Earth's bow shock, models of supernova remnant 
shocks, and the heliospheric termination shock (e.g., Ellison, 
Jones & Baring 1999). Additionally, while this paper worked 
primarily with large Alfvenic Mach number shocks, in princi- 
ple, low Alfvenic Mach numbers are possible in jet systems as 
well. While the Rankine-Hugoniot solutions presented above 
are fully capable of determining the appropriate jump condi- 
tions, low Alfvenic Mach number shocks may produce signifi- 
cant second-order Fermi acceleration due to the motion of the 
scattering centers (Alfven waves) in the upstream and down- 
stream rest frames. The code currently assumes scattering cen- 
ters that are stationary in their respective fiuid frames but can 
easily be adapted to include non-stationary scattering centers 
for the case of low Alfvenic Mach number shocks. In addition, 
preliminary work has also been done laying the ground work for 
future inclusion of cross-shock potentials in the simulation (Bar- 
ing & Summerlin, 2007). The simulation is currently a single 
fluid model, treating electrons/pairs or ions. For electron-proton 
shocks, the disparate diffusion scales of the two species will cause 
their distribution functions to react to the presence of the shock 
on different length scales. This charge separation at the shock 
discontinuity induces an electric field that acts to restore quasi- 
neutrality, and can lead to significant energy exchange between 
ions and electrons. The inclusion of these effects and the deter- 
mination of their impact on injection and acceleration of protons 
and electrons will be the focus of future work. 
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APPENDIX 



APPENDIX A: THE ANGULAR DISTRIBUTION OF PARTICLES AT THE RETURN PLANE 

Downstream of the shock, once the particle distribution has realized isotropy in the local fluid frame at some position x , it 
maintains isotropy as an asymptotic state at all positions further downstream. Tracking diffusion downstream of x to assess the 
momentum components of particles that return to position x is CPU-intensive. A much faster method is to compute the probability 
of return to x and these momentum components statistically, subject to the condition of isotropy in the local fluid frame, and 
constancy of the returning particle's fluid frame momentum pf . This was the expedient approach of Ellison et al. (1990a). In 
this Appendix we develop the formalism for angular distributions of returning particles for arbitrary pf , not just ultra- relativistic 
particles, as has been the restriction of previous expositions. At x , the flow speed will be /3c , i.e. of Lorentz factor F (subscripted 
2 in the main text). The particle will have a dimensionless momentum P/=7//3/ (Ps=7s/3s) and angle cosine with respect to 
the shock normal of /i^ (fig) in the fluid (NIF) frame. The non-covariant momentum distribution function downstream of the 
probability of return plane located at x (and all positions further downstream) assumes the form 



NSjpf -po) 



Anp 



f 



(Al) 



in the local fluid frame. Observe that N , the total number of particles, is a Lorentz invariant under boosts along the shock 
normal, and / is to be distinguished from the phase space density /(r, p) . Equation 10.2 of Landau and Lifshitz (1989) then 
gives the transformation of the distribution functions as "fsfsiPs) — IfffiPf) ! in covariant formulations, the particle Lorentz factor 
is absorbed into the definition of / . This is employed in the second integral in Eq. (|Al[) , along with the fluid frame distribution 
function. To integrate over the delta function, it is necessary to change variables dps dpf , so that 
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Hereafter, the identity Pf = po will be assumed. The partial Jacobian \dps/dpf\ can be determined by first writing the Lorentz 
boost relations for the particle momentum: 
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After moderate algebraic manipulation, one can eliminate and solve for Ps/pf as a function of pf and '■ 
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Additionally, holding the integration variable, pLg , constant, one can derive 
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These allow us to rewrite Eq. (|A2p in the form 
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The angular distribution dN/dfis describes the beaming appearing in the NIF shock frame of an isotropic distribution in the 
fluid frame of fixed, specified momentum pj . It is applicable to arbitrary 7/ , not just ultra-relativistic cases 7/^1, the usual 
restricted consideration (e.g. see Peacock 1981). The identity for dN/dfis , casting it in terms of a perfect derivative dJ^/dfi^ , can 
be established using a moderate amount of algebra; it nicely facilitates the derivation of the integral identity for N in Eq. (jA6p . 

To formulate probabilities of transmission and return at position x , the angular distribution dN/dfig must be weighted by the 
flux of particles incident upon the plane at x that is parallel to the shock plane. In the NIF, this weighting factor is proportional 
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to the density of partieles, which is proportional to 7^ , and also to the velocity component Ws|/^s| along the a; -direction. In this 
way, we form flux angular distributions using Eq. (|A4p as follows: 



= c\^is\4P^^l#t^ ^ , (A7) 

dus ' S'(l -^Vs) ' ' 

for C being a constant of normalization. This formula generalizes that employed in Peacock (1981), which is restricted t o 7/ ^ 1 
cases. It is easily seen that the integrands in Eq. (26) of Peacock's paper are proportional to the 7/^-00 limit of Eq. (|A7|) . and 
can essentially be derived using light aberration considerations. The total probability Pr of return to x of particles of fixed pf 
incident from the upstream side is then the ratio of two integrals over the flux distribution: 

This simple expression for Pr is valid for any flow speed /3 and any particle speed /?/ > ; it was flrst derived for non-relativistic 
shocks (i.e. /3 ^ 1 ) by Bell (1978) and for relativistic shocks with high-speed (/?/ « 1) particles by Peacock (1981). Again, 
recognizing the appearance of a perfect derivative d{Y?) j d^s in the flux distributions expedites the integrations in Eq. (jASp . The 
direction of the magnetic field, and the level of cross-field diffusion are irrelevant to the derivation of the flux angular distribution 
and probability of return. Technically these formulae, as derived, must be applied in a downstream normal incidence frame] in 
practice, for many Rankine-Hugoniot solutions such as for high Alfvenic Mach numbers (e.g., see Figs. [2] and [3]), the flow deflection 
at the shock is small, and the downstream and upstream NIF frames are almost coincident. 

The simulation computes the statistical probability of return to the plane at x according to Pr ■ Those particles that are 
deemed to escape are eliminated from the simulation. For those that return, their returning value of < is selected randomly 
from the distribution in Eq. (jA7[) when it is normalized to unity on — 1 < /^s < 0. Then the constant of proportionality is 
C = 2jjP'j/r^/{j3f — /3)^ . Upon return, the particle is placed at x with the same (y, z) coordinates it originally crossed with; the 
system is uniform in the dimensions transverse to the shock normal, so this step introduces no bias. After /^s is determined, the 
phase of the momentum vector about the shock normal is selected randomly from a uniform distribution on [0, 27r] . The returning 
particle momentum vector is then totally specified, and is routinely cast in rotated coordinates to identify variables connected to 
gyration about oblique magnetic fields. 



